{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Variacional Univariate Gaussian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from numpy.random import randn, normal, seed\n",
    "from scipy.stats import norm, gamma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "%config InlineBackend.figure_format = \"retina\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook we look at a factorized variational approximation to a Gaussian Distribution.\n",
    "\n",
    "Let $\\mathcal D= \\{x_n\\}_{n=1}^N$ i.i.d. random variables such that $x_n\\sim \\mathcal N(\\mu, \\tau^{-1})$. Our goal is to find a posterior distribution \n",
    "\n",
    "$$\n",
    "    p(\\mu, \\tau | \\mathcal D) \\propto p(\\tau)p(\\mu|\\tau) p(\\mathcal D | \\mu, \\tau)\n",
    "$$\n",
    "\n",
    "With\n",
    "* $p(\\tau) = \\text{Gam}(\\tau|a_0, b_0)$\n",
    "* $p(\\mu|\\tau) = \\mathcal N(\\mu|\\mu_0, (\\lambda_0\\tau)^{-1})$\n",
    "* $p(\\mathcal D | \\mu, \\tau) = \\prod_{n=1}^N \\mathcal{N}(x_n | \\mu, \\tau^{-1})$\n",
    "\n",
    "To find such posterior, we assume a factorized variational approximation to $p$ given by\n",
    "\n",
    "$$\n",
    "    q(\\mu,\\tau) = q^*_\\mu(\\mu) q^*_\\tau(\\tau)\n",
    "$$\n",
    "\n",
    "\n",
    "Since the optimum factors of a factorized variational approximation follow the relation\n",
    "\n",
    "$$\n",
    "    \\log q^*_j({\\bf z}_j) = \\mathbb{E}_{m\\neq j}[\\log p({\\bf z})] + C,\n",
    "$$\n",
    "\n",
    "It follows that\n",
    "\n",
    "$$\n",
    "    q^*_\\mu(\\mu) = \\mathcal{N}\\left(\\mu \\ \\big | \\frac{N\\bar x + \\mu_0}{N + \\lambda_0}, [N + \\lambda_0]\\mathbb{E}[\\tau]\\right)\n",
    "$$\n",
    "\n",
    "\n",
    "$$\n",
    "    q^*_\\tau(\\tau) = \\text{Gam}\\left(\\tau \\ \\big | a_0 + \\frac{N}{2}, b_0 + \\frac{1}{2}\\mathbb{E}\\left[\\lambda_0\\tau(\\mu - \\mu_0)^2 + \\sum_{n=1}^N(x_n - \\mu)^2\\right] \\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x1a26f84490>"
      ]
     },
     "execution_count": 147,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAvUAAAHwCAYAAAAmZ5CjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAWJQAAFiUBSVIk8AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeXjcZ333+893tO+LJdmyZNmWd8e7nTj7TghJSCiQ9jRA01Da8kAPPDzwtL3KdcoDp5zTXi2UwgOlpy0EytImAUJIAgSykcQhifd9t2zJkiVLsvZ95j5/zGhmpFi2ZI30m+X9ui5fM7+Z38zvKy6i+eie733f5pwTAAAAgMTl87oAAAAAANNDqAcAAAASHKEeAAAASHCEegAAACDBEeoBAACABEeoBwAAABIcoR4AAABIcIR6AAAAIMER6gEAAIAER6gHAAAAEhyhHgAAAEhwhHoAAAAgwaV7XUC8M7NTkgol1XlcCgAAAJLbIkldzrnFU30hof7yCnNyckpXrVpV6nUhAAAASF6HDh1Sf3//Fb2WUH95datWrSrdsWOH13UAAAAgiW3evFk7d+6su5LX0lMPAAAAJDhCPQAAAJDgCPUAAABAgiPUAwAAAAmOUA8AAAAkOEI9AAAAkOAI9QAAAECCI9QDAAAACY5QDwAAACQ4Qj0AAACQ4Aj1AAAAQIIj1AMAAAAJjlAPAAAAJDhCPQAAAJDgCPUAAABAgiPUAwAAAAmOUA8AAAAkOEI9AAAAkODSvS4AAKbqB2+c8bqEuPXQ1hqvSwAAeICRegAAACDBEeoBAACABEeoBwAAABIcoR4AAABIcDEL9WZWbWbfMrNGMxs0szoz+4qZlUzhPf6nmT0bem2PmXWZ2T4z+7KZVV/idavN7DEzazGzATM7YmafN7Oc2Px0AAAAQPyKyeo3ZrZE0jZJFZJ+KumwpGskfVLS3WZ2g3OubRJv9aeSeiS9LKlZUoakjZI+JemPzOxW59yucdfeKumF0LlPSKqXdLukv5Z0h5nd4ZwbnP5PCQAAAMSnWC1p+Q0FA/0nnHNfG33QzL6sYCD/oqSPTuJ91jjnBsY/aGZ/LOn/C73PPVGPp0n6tqRcSQ84554KPe6T9Jik94Wu/7dX9mMBAAAA8W/a7TdmVivpLkl1kr4+7unPSeqV9CEzy7vce10s0Ic8FrpdNu7xWyStkvSb0UAfep+ApD8PHX7UzOxy1wYAAAASVSx66m8P3T4XCtNhzrluSa8pOJJ+7TSu8e7Q7d4Jrv2L8S9wzp2UdFTSQkm107g2AAAAENdi0X6zInR7dILnjyk4kr9c0vOTeUMz+4ikakn5ktZKulPSaUl/eQXXXh76d+Iy19wxwVMrJ1EyAAAA4JlYhPqi0G3nBM+PPl48hff8iKStUcdvSXrIOXd8Fq4NAAAAJJRYTZS9lNF+djfZFzjnrpUkM5sjaZOCE2R3mNnvOefe1moTi2s75zZf9A2CI/ibpnBNAAAAYFbFoqd+dDS8aILnC8edN2nOuTbn3K8UbN/pl/TdcWvPz9i1AQAAgEQRi1B/JHS7fILnR1esmajv/bKccx2SXpdULumq2bw2AAAAEO9iEepfDN3eFVofPszMCiTdoOAo+2+neZ2q0O1I1GMvhG7vHn9yaKnN5QpOsD05zWsDAAAAcWvaod45d0LSc5IWSfr4uKc/LylP0nedc72jD5rZSjMbs6qMmS0MBfG3MbM/lXS1grvF7ot66mVJhyTdbGb3R53vk/R3ocNvOucm3c8PAAAAJJpYTZT9mKRtkr5qZncoGLS3SrpNwdaXz447/1DoNnpTqI2Sfmxm20KvaZY0R8H17ddK6pH0Ieecf/QFzjm/mT2i4Ij9E2b2hKQzku6QtEXBNfL/MUY/IwAAABCXYtF+Mzpav0XSowqG+U9LWiLpq5Kuc861TeJtdioYwDMl3SvpM5J+X8GVa74kabVz7uWLXPsNBUfxf6rghNpPKThx9guS3uGcG5zOzwYAAADEu5gtaemcq5f0yCTPtYs8dkbBPwau5NoHJT14Ja8FAAAAEl1MRuoBAAAAeIdQDwAAACQ4Qj0AAACQ4Aj1AAAAQIIj1AMAAAAJjlAPAAAAJDhCPQAAAJDgCPUAAABAgiPUAwAAAAmOUA8AAAAkOEI9AAAAkOAI9QAAAECCI9QDAAAACY5QDwAAACQ4Qj0AAACQ4Aj1AAAAQIIj1AMAAAAJjlAPAAAAJDhCPQAAAJDgCPUAAABAgiPUAwAAAAmOUA8AAAAkuHSvCwAAxM4P3jjjdQlx6aGtNV6XAAAzipF6AAAAIMExUg8AuKSAc+rqH1b/sF/OSc4FH3POKTsjTcW5mcpMZ4wIALxEqAcASJJG/AE1dQ7oTHufznUN6ELfkDr6htXRN6SAu/RrczPTVJKbqeLcDFUV52hpRb7mF+fIZzY7xQNAiiPUA0CKGvEHdLK1V8eau1V/oV9nO/rlv1x6n0DfkF99Q8H3ONDYpecONisnI01LK/K1tCJfV80vVG4mHzkAMFP4DQsAKaR3cERHznXr0LkuHWvp0dBIYFKvy8tMU15WutJ8JpNkZjILvl9n//BFR/L7h/3ad7ZT+8526um9jdpYU6Lra+eoojA7tj8UAIBQDwDJbtgf0OFz3dp5+oKONnfrUmPxpXmZqinNVXVJjubkZakkN+OyPfMB59Q9MKKOviG19gzqxPleHW/pUc/gSFQNTm+eatebp9q1fG6+rl9SpmUV+TLacwAgJgj1AJCkzl7o1/bT7drb0Kn+Yf9FzynNy9SqeQWqLc/XgtJc5WdN/WPBZ6ainAwV5WRo4Zw8bV5YKuecznUN6HhLj3bXd6ipcyB8/tHmHh1t7tGS8jw9sL5KZQVZV/wzAgCCCPUAkET8AaeDTV167XirzrT3XfScmtJcra4s1Mp5BSovyJqR0XIzU2VRjiqLcnTj0jKdau3VayfadLipK/xNwYnzvfqnF47p1uXlumV5udLTWEEHAK4UoR4AksDAsF/b69q17WSbOvqG3/Z8SW6GNtWUaGNNiUrzMme1NjNTbXm+asvz1dYzqFePt+rNU+1yCv4R8vzhFu1p6NADG6q0pDx/VmsDgGRBqAeABDYw7NdrJ1r12vFWDQyPnfSaZqa11UXasrBEi8ry4mJ5yTn5WXpgQ5WuXlSqJ3efVcOFfklSa8+Q/v3VU7p1RbnuXDU3LmoFgERCqAeABDQ44tfrJ9r0yrHWt/XL52amaeviUm2tnaPC7AyPKry0+cU5+ugtS/TGqXY9d+CcBkOr8Lx05Lyauwb1u5urlZWR5nGVAJA4CPUAkEBGAgG9cbJdLx5pUd/Q2DA/Jy9TNy0r18aaYmUkQH+6z0zX1c7RVZWF+tHOBh1r6ZEkHWrq0j+/fEIfunah5uQziRYAJoNQDwAJ4si5bj2zr0mtPYNjHi/Ny9TtKyu0vrpYab7Ea1spzMnQw9cv0i/3n9Mrx1slSS3dg/rGSyf00NYa+uwBYBII9QAQ51p7BvXM3iYdae4e83hxboZuX1GhjTUlCRnmo/nM9K61lZpblK2f7Dorf8Cpf9ivb792Sr9/TY2uml/kdYkAENcI9QAQp0b8Ab14pEW/Odoqv4tsGZWV7tPtKyt03ZI5SvfFf5vNVGyqKVFZfpa+/9vT6h4cUcBJ//lmvT6w1bSystDr8gAgbiXXpwEAJImGC336+kvH9eKR8+FAb5K2LCzR/3jHct20rDzpAv2omtJcfey2pZoTWnrT75y+/+YZHR33TQUAICI5PxEAIEGN+AN67sA5ffPlE2ruivTO15Tm6mO3LtV7N1WrIE5XtImlopwMfeSmWpXkBn9Wf8Dpe789rRPnezyuDADiE6EeAOJEU2e//veLx/XS0fMKhLptMtJM715XqT+5uVZVJTneFjjLinIy9JEba1WcEwz2IwGn775ep1Otvd4WBgBxiFAPAHFge127/vmlE2rpjozOLy7L0yduX6brlpSl7GZMJXmZ+qMbF6swOzgFbNgfDPYt3QPeFgYAcYZQDwAeGhoJ6IkdDfrxrrMaCQ3Pj47O/9GNi1mnXcFdaP/oxloVZAWD/eBIQN/77Wn1j1unHwBSGaEeADzS2j2ob758QjvPXAg/NrcwS392W2qPzl9MeUGWHr5+kTLSgv+btPYM6bHt9QpErQoEAKmMUA8AHjja3K2vv3Rc57oibSQbFxTrv92yVOUFjM5fzPziHL1vU3X4+Ehzt359sNnDigAgfrBOPQDMsu117Xpy99nwZNh0n+nd6+Zry6ISGaPzl7SuuliNHQP6zbHzkqSXjp5XZXGO1laxORWA1EaoB4BZ4pzTC4db9PzhlvBjxTkZ+uC1CzW/OLVWtpmOu66aq3Nd/TraHFze8okd9SrLz1RlEf8bAkhdtN8AwCzwB5x+suvsmEBfWZStj966hEA/RT4z/d6WmvDmVMN+px+8cUZDIwGPKwMA7xDqAWCGDfuDq7VsPx2ZELu0Il9/fFOtClNgI6mZkJOZpg9eu1CZ6cGPsbbeIf18f5PHVQGAdwj1ADCDRvwBff+N0zrS3B1+bFNNsR6+bpGyM9I8rCzxzS3M1v3r5oeP3zjVrmNR/zsDQCoh1APADBkJBPSDN8+Ee78l6dYV5Xrfpmql+ZgQGwsba4q1urIwfPyjnQ2sXw8gJRHqAWAG+ANO//VWvQ6fi4wc376yQnetnscKNzFkZnrPxirlZQa/9egaGNHP9jZ6XBUAzL6YhXozqzazb5lZo5kNmlmdmX3FzEom+fo8M/uAmf3AzA6bWa+ZdZvZdjP7tJllTvA6d4l/v43VzwcAkxVwTo/vqNeBxq7wY7csL9cdKys8rCp55Wel63c2Rtav313foX1nOz2sCABmX0yWtDSzJZK2SaqQ9FNJhyVdI+mTku42sxucc22XeZubJH1PUrukFyU9KalU0rsl/YOk95rZHc65gYu89rSkRy/yeMPUfxoAuHLOOf14Z4P2NkRC5Y1Ly3TX6rmM0M+g1fMLtammJLw77093n9XCOblMRAaQMmK1Tv03FAz0n3DOfW30QTP7sqRPSfqipI9e5j3OSfqgpMedc0NR71Eg6SVJ10v6uKQvXeS1dc65/zWN+gEgJl480qKdZzrCx9fWztG71tByMxvuW1epk+d71NE/rL4hv362p1Ef2LrQ67IAYFZMu/3GzGol3SWpTtLXxz39OUm9kj5kZnmXeh/n3G7n3PejA33o8W5Fgvyt060XAGbK3oYO/fpQZB36LQtLdN+6SgL9LMnOSNP7NkfacA40dukoq+EASBGx6Km/PXT7nHNuzM4foUD+mqRcSddO4xrDoduRCZ4vNrMPm9lfmdnHzWw61wKAKWu40KcndkQ6/paU5+mBDVXyEehn1ZLyfG2qKQ4f/2xPo0b8bEoFIPnFov1mRej26ATPH1NwJH+5pOev8BofDt3+YoLn10v69+gHzGyPpA855/ZN5gJmtmOCp1ZOqkIAKauzf1j/8dvTGgk4SVJZfqYeumYhy1Z65J1XzdPBpi4NDAfU1jukV4+36g+uX+R1WQAwo2IxUl8Uup1oqYHRx4sneP6SzOzPJN0tabekb13klC9LukFSuaQCSVdLekLBoP+CmVVdyXUBYDKGRgL6j9fr1D0Q/CIxJyNNf3DdIuVksrGUVwqyM3Tnqrnh4xePtOhsR7+HFQHAzJuNdepHh6rclF9o9l5JX1FwEu37nHPD489xzn3aObfNOdfqnOtxzm13zj0o6UeSyiR9ZjLXcs5tvtg/BVfyAYC3cc7piZ0NauwMLsrlM+mhrTUqy8/yuDJsXTxH8wqzJUnDfqcvPnPQ44oAYGbFItSPjsQXTfB84bjzJsXM3iPpPyW1SLrVOXdyinV9M3R78xRfBwCT8lbdBe2PWg/9gfVVWlKe72FFGJXmM92/fn74+Nl95/TqsVYPKwKAmRWLUH8kdLt8gueXhW4n6rl/GzN7UNLjkpol3eKcO3KZl1zM+dDtJVfdAYAr0dw1oGf2RXYuvWZxqa5eXOphRRhvUVmeNi6IdH7+9VP7NTTCpFkAySkWof7F0O1dZjbm/UJrzN8gqV/SpHZ3NbOHJP1QUqOCgf7YFdY1ugLOVEf4AeCShv0B/ddb9Rr2B7sK5xZm6d61lR5XhYu5e808ZaUHP5pOnu/Vd7bVeVsQAMyQaa9+45w7YWbPKbjCzcclfS3q6c8rOFL+L8653tEHzWxl6LVj+tXN7GEFJ8OelnSbc+70pa5tZpskHYl+79Dj6xTc8EoK7lILADHzi/3ndK4r2Eef7jP93tU1ykibjSlKmKrRSbPP7GuSJH35V0eV5jNlZzCRedRDW2u8LgFADMRqR9mPSdom6atmdoekQ5K2SrpNwbabz447/1DoNrzem5ndpmCg9yk4+v/IRTZs6XDOfSXq+BOS3mtmL0iqlzSo4BKUd0tKk/SvCo76A0BMHGrq0usn28LH96ytDE/IRHzaWluq10+2qb13SP3Dfr1y7LzesXqe12UBQEzFJNSHRuu3SPqCgoH6HklNkr4q6fPOufZJvM1CRdqBPjzBOacVXA1n1JMKTsRdp+AmWNmS2iT9XNK/OueemuKPAgAT6uwf1o92RjaYWl1ZqK300ce9dJ9Pd66aq8e210uSXj3eqmtr56ggO8PjygAgdmI1Ui/nXL2kRyZ57tuG4J1zj0p6dIrXfFLBYA8AM8o5p5/salDfkF+SVJidrvdurNJFvlFEHFpXXaTfHD2vc10DGvY7vXikRfevZxsTAMmDJlAAmIS9Zzt1tLlHUrBv8He3LFBuVszGRTDDfGZ651WRDanePNWu9t4hDysCgNgi1APAZfQP+fX03qbw8dbaOaplPfqEs3xugRbNyZUkBZz060PNHlcEALFDqAeAy/j5/ib1Do5ICrbd3LV67mVegXhkZnrnVZEJsnvqO9TU2e9hRQAQO4R6ALiEk6092n76Qvj4/vXzWQ4xgS2ck6eV8wokSU7ScwcYrQeQHAj1ADCBEX9AT+6K7Bq7urJQq+cXeVgRYuGu1fPC6ykfae7W6bbeS54PAImAUA8AE3jp6Hm19gxKkrLSfXr3+vkeV4RYmFeUrQ0LisPHLx5p8bAaAIgNQj0AXERL14BePnI+fPzOq+apKId1zZPFrSsqwqP1R5t7dLaD3noAiY1QDwAX8ez+JvmdkyTVlObqGjaZSirlBVlaUxVppXqJ0XoACY5QDwDjHGvpHrMm/QMb5svHJlNJ59YV5eH7Bxq71Nw14GE1ADA9hHoAiBJwTj/fdy58vHlhiSqLcjysCDOlsignvBKOJL189PwlzgaA+EaoB4Aou85c0LnQiG1GmunOVaxJn8xuXVERvr+3oYNdZgEkLEI9AIQMjQT03MHIuuU3LytXIZNjk1pNaa5qy/MkBXeZZbQeQKIi1ANAyCvHz6t7ILhzbEF2um5aVn6ZVyAZ3BY1Wr/zzAV19g97WA0AXBlCPQBI6hoY1itHW8PH71g1V5np/IpMBbVleaopzZUk+QNOrxxjtB5A4uETCwAk/fpgs4b8AUnSvMJsbVpY4nFFmC1mNmYlnLfq2tU3OOJhRQAwdYR6ACmvuWtAO05fCB+/a808lrBMMSvmFqiyKFuSNOx3erOu3eOKAGBqCPUAUt4Lh1vkQveXVeRr2dyCS56P5GNmumFpWfj49ZNtGgkEPKwIAKaGUA8gpTV3DWj/2c7w8TtWs4RlqlpXXaSCrHRJUvfAiPY1dF7mFQAQPwj1AFLai0cio/Qr5haouiTX03rgnXSfT9cumRM+fu1Eq5xzl3gFAMQPQj2AlNXSPTBmNPb2lRWXOBup4JpFpUr3BedTNHYMqK6tz+OKAGByCPUAUtZLR86HR+mXz83XglJG6VNdXla6NtZEVj567XjrJc4GgPhBqAeQklq7B7WnviN8fPtKeukRdENUC86hpi619Qx6WA0ATA6hHkBKiu6lX1aRH958CKgozNbyufmSJCdp28k2bwsCgEkg1ANIOW09g9rTED1KTy89xrphSWR5yx2nL2hg2O9hNQBweYR6ACnnxSPnFQgN0y8tz9fCOXneFoS4s7QiXxUFWZKkoZGA3mIzKgBxjlAPIKV09A1pd31k91hG6XExF9uMKsDylgDiGKEeQErZdqItPEq/uCxPi8oYpcfFbVhQrNzMNElSR9+wjjZ3e1wRAEyMUA8gZQwM+8e0Udy8rOwSZyPVZaT5tHlhZHnLN07SggMgfhHqAaSM7XXtGhwJSJLK87O0bG6BxxUh3l2zqDR8/2hzt9p7hzysBgAmRqgHkBL8AadtJyJLE964tEw+Mw8rQiKYk581ZnlLJswCiFeEegApYX9jpzr6hyVJeZlp2lBT7HFFSBRbF0c2o3qrrl0j/oCH1QDAxRHqASQ955xePdYaPt5aO0cZafz6w+SsmFegopwMSVLfkF/7G7s8rggA3o5PNQBJr66tT2c7+iVJ6T7TtbVzLvMKIMJnpqujeuvfOMUOswDiD6EeQNJ79XhklH5jTbHys9I9rAaJ6OpFJfKFpmCcbuvTuc4BbwsCgHEI9QCSWmvPoA43RdolbljCMpaYuoLsDF01vyh8zGg9gHhDqAeQ1F473qrRfUBXzC1QRWG2p/UgcW1dHGnB2VXfocFhv4fVAMBYhHoASWtg2K+dZy6Ej29ksylMw+KyPJUXZEmShkYC2lXf4XFFABBBqAeQtHaduaBhf3Ccfm5hlmrL8jyuCInMzMaM1m9nzXoAcYRQDyApOef021OR0LV18RwZm01hmjYuKFF6aMZsY+dAeFUlAPAaoR5AUjrV1qvz3YOSpMw0nzYsYLMpTF9OZprWVEUmzDJaDyBeEOoBJKU3TkbC1oYFxcrOSPOwGiSTLQtLwvf3NHRomB1mAcQBQj2ApNM9MKyDUbt+bq0tvcTZwNQsLstTaV6mJGlgOKADjZ0eVwQAhHoASWjH6Qvyu+AE2ZrSXFUW5XhcEZKJmY0Zrd9ed+ESZwPA7CDUA0gqAef05pgJsozSI/Y21pRodNr1ydZetfUMeloPABDqASSVI+e61dE/LEnKHTepEYiVopwMrZhXED7ecZrRegDeItQDSCpvnGoL39+ysEQZafyaw8yIbsHZeeaC/AF3ibMBYGbxaQcgabT3DulYc48kySRds3iOtwUhqa2YV6j8rHRJUtfAiI61dHtcEYBURqgHkDTePNWu0bHSZXPzwyuUADMhzWfaWBPZ/4AJswC8RKgHkBT8AaedZyKhaiuj9JgFWxZGJmIfPtel7oFhD6sBkMoI9QCSwtHmbvUMjkiSCrLTtXxuwWVeAUxfeUGWFs7JlSQFnLS7vsPjigCkKkI9gKQQvfrIppoSpfnsEmcDsTN+wqxzTJgFMPsI9QASXvfAsA6fi+wgu7mm5BJnA7G1Zn6RMtKCf0Q2dw2qsXPA44oApCJCPYCEt7u+Q6OrCS6ck6uygixvC0JKycpI01XzI/sh7DrDhFkAsy9mod7Mqs3sW2bWaGaDZlZnZl8xs0kNmZlZnpl9wMx+YGaHzazXzLrNbLuZfdrMJlzGwsxWm9ljZtZiZgNmdsTMPm9m7A0PJDnn3JjWm+hWCGC2bIr6dmhPfQdr1gOYdTEJ9Wa2RNIOSY9IelPSP0o6KemTkl43s8ksQ3GTpO9Jeqek/ZK+JumHkqok/YOkF80s+yLX3irpLUnvkfRrSf8kqUvSX0v6lZkxZAcksYYL/WrpHpQkZab52EEWnqgtz1NhdnDN+t4hv442s2Y9gNkVq5H6b0iqkPQJ59x7nHN/6Zy7XcFwv0LSFyfxHuckfVBSpXPu/aH3+BNJyyXtlHS9pI9Hv8DM0iR9W1KupPc75x5yzv2FpK2SfiTpBkmfislPCCAubY8apV9bXaSs9DQPq0Gq8plpY83YCbMAMJumHerNrFbSXZLqJH193NOfk9Qr6UNmlnep93HO7XbOfd85NzTu8W5JXwod3jruZbdIWiXpN865p6JeE5D056HDj5oZy2AASWhoJKC9DZElBGm9gZc2LohsRHW4qVt9QyMeVgMg1cRipP720O1zoTAdFgrkryk4kn7tNK4xupvH+N+Qo9f+xfgXOOdOSjoqaaGk2stdwMx2XOyfpJXTqBvADDrQ2KnBkeCvnbL8TNWU5npcEVJZRWG2qkuCU7n8zmlvQ6fHFQFIJbEI9StCt0cneP5Y6Hb5NK7x4dDt+PA+G9cGEKeiW28215SIL+XgNVpwAHglPQbvMTorbaIhidHHiyd4/pLM7M8k3S1pt6RvzdS1nXObJ7j+DkmbLl8pgNnU1jOoU629kiTT2DAFeGV9VZGe3dskv3OhSdwDqih42xoPABBzs7FO/ejQ2ZTX9zKz90r6ioKTaN/nnBu+zEtidm0A8W3nmUgv/fK5BSrMyfCwGiAoNytdK+YVhI93Rf3/FABmUixC/eho+ETryBWOO29SzOw9kv5TUoukW0M98rNybQDxzTmn3fWR1oZNTJBFHIlesz64MRrjSgBmXixC/ZHQ7UR968tCtxP1vb+NmT0o6XFJzZJucc4dmeDUmF8bQPw7096nC33BL+6yM3xaGTUyCnht+bx85WYGl1bt7B/WyfO9HlcEIBXEItS/GLq9y8zGvJ+ZFSi4Vny/pN9O5s3M7CEFN51qVDDQH7vE6S+Ebu++yPvUKhj2Tyu4ERaAJBHd0rC2qkgZabPRSQhMTrrPp/XVkalce+ppwQEw86b9SeicOyHpOUmLNG5zKEmfl5Qn6bvOufBQhZmtNLO3LRVpZg9L+g9JZyTdPEHLTbSXJR2SdLOZ3R/1Pj5Jfxc6/KZzfPcJJIvBEb/2nY101G1YQOsN4s+GqDXr9zd2atgfuMTZADB9sVj9RpI+JmmbpK+a2R0KBu2tkm5TsPXls+POPxS6Da8/Z2a3Kbi6jU/B0f9HLrI8XYdz7iujB845v5k9ouCI/RNm9oSCfxDcIWmLgmvk/2MsfkAA8eHFw+fVP+yXJBXnZmjhHNamR/ypLsnRnLxMtfUOaXAkoMPnurW2aqLpXwAwfTEJ9c65E2a2RdIXFGyFuUdSk6SvSvq8c659Em+zUJFvDj48wTmnFVwNJ/rab5jZ1Qp+K6+nUZIAACAASURBVHCXpILQeV+Q9LfOucEp/jgA4thPdjWE729YUCwfa9MjDpmZ1i8o1guHWyRJu89cINQDmFGxGqmXc65e0iOTPPdtn8LOuUclPXqF1z4o6cEreS2AxNHRN6QXD58PH0e3OADxZkNUqD/a3KO+wRHlZsXsYxcAxmB2GYCE8cy+Jg2FepOrinPY1AdxrSw/S9UlOZIkv3Pa18jqygBmDqEeQML4yc6z4fsbaxilR/yL/jZpN6vgAJhBhHoACeFMW5+2nw5uOOUzaV01oR7xb111sXyhhtPTbX260DvkbUEAkhahHkBCeHJ3ZJR+WUWB8ulNRgLIz0rX0or88PGeBkbrAcwMQj2AuOec0092RUL9BlpvkECiW3B21XeIrVMAzARCPYC4t6ehU6dag/vX5Wela9W8Qo8rAiZvVWWhMtKCPTjnuwfV1DngcUUAkhGhHkDce2p3Y/j+O6+ap8x0fnUhcWSlp+mq+ZE16pkwC2Am8MkIIK75A05P742E+gc2zPewGuDKrI+a2L2noUMBWnAAxBihHkBce+NUm1q6gxtDl+Vn6volczyuCJi6pRX5ygtN7u4eGAm3kwFArBDqAcS1n+2JjNLfs7ZS6Wn82kLiSfOZ1lZF5oLsZRUcADHGpyOAuDU0EtCz+86Fj+9fT+sNEld0C87+s10aCQQ8rAZAsiHUA4hbrxw7r87+YUlSVXGONtWUeFwRcOUWlOaqKCdDktQ/7Nfxlh6PKwKQTAj1AOLWU1GtN/etr5RvdGtOIAH5zLSuOrIKzt6GTg+rAZBsCPUA4lL/kF+/OtgcPqb1BslgXVQLzsGmLg2N0IIDIDYI9QDi0q8PNatvyC9JWlKep9WVbDiFxDe/KFtl+ZmSgnNGjjR3e1wRgGRBqAcQl6Jbb+5fXyUzWm+Q+MxszGg9q+AAiBVCPYC409k/rJePnA8f38+GU0gi0X31R851a2DY72E1AJIFoR5A3Pnl/nMa8gd7jddWFWlxWZ7HFQGxU1GQrcqibEnSSMDpYFOXxxUBSAaEegBxJ7r15t3rKz2sBJgZtOAAiDVCPYC4cr57UNtOtIaP71tH6w2Sz7qqSAvO8ZYe9Q6OeFgNgGRAqAcQV35x4JwCLnj/6kUlml+c421BwAwoyctUTWmuJCngpP2NrFkPYHoI9QDiyjN7ozacYpQeSYyNqADEEqEeQNxo6R7QG6faJUlm0rvWzPO4ImDmrK0q0uhCrXWtversH/a0HgCJjVAPIG78Yv85uVDrzTWLSlVRmO1tQcAMKsjOUG15cGUnJ2kfE2YBTAOhHkDceHpvU/j+fetY9QbJb330KjhnacEBcOUI9QDiQnPXgN6qC7be+Ey6ew2hHsnvqvlFSgvtltxwoV9tPYMeVwQgURHqAcSFn+9rCrfeXFs7R+UFWd4WBMyCnMw0LZubHz5mtB7AlSLUA4gLz+yLtN7cS+sNUkh0C86eevrqAVwZQj0Az53rHNBbdRckhVpvrmLVG6SOlZUFykgLtuC0dA/qXNeAxxUBSESEegCeix6lv35Jmebk03qD1JGVnqaV8wrDx3sZrQdwBQj1ADwXveEUrTdIReujN6I62yk3OsEEACaJUA/AU2c7+rXzTHBkMs1neietN0hBy+cWKDsj+JHc3jukhgv9HlcEINEQ6gF46udjWm/mqDQv08NqAG+kp/m0ujJqtJ6NqABMEaEegKfYcAoIGt+CE6AFB8AUEOoBeKa+vU+7Q5MC032mu1bTeoPUVVuer7zMNElS98CI6lp7Pa4IQCIh1APwzM/3R0bpb1haphJab5DC0nymNVWR0fo9DWxEBWDyCPUAPPMMrTfAGNEbUR1o7JQ/QAsOgMkh1APwRH17X3gkMiON1htAkmrm5KooJ0OS1Dfk14nzPR5XBCBREOoBeCJ6w6mblpWrKDfDw2qA+OAz09oqVsEBMHWEegCeiG69uXctrTfAqHVRq+AcaOzSiD/gYTUAEgWhHsCsO93Wq31ng603mWk+3bl6rscVAfGjqjgnvF/D4EhAR5tpwQFweYR6ALMuem36m5eXhXuIAUhmpnXRLThnacEBcHmEegCzbkzrDaveAG+zLmoVnENNXRoaoQUHwKUR6gHMqpPne3SwqUuSlJnu052raL0BxptbmKWKgixJ0rDf6fC5Lo8rAhDvCPUAZtWzUave3LK8XAXZtN4A45nZmAmze9mICsBlEOoBzKqn2XAKmJR1VZEWnCPN3RoY9ntYDYB4R6gHMGuOt/To8LluSVJWuk930HoDTKisIEvzi7IlSf6AC7etAcDFEOoBzJro1pvbVlQoPyvdw2qA+Bc9YZaNqABcCqEewKyJXvXmHlpvgMuK3l32eEuPegdHPKwGQDwj1AOYFceau3WkOdh6k53h0x0rKzyuCIh/JXmZqinNlSQFXHCHWQC4GEI9gFkRPUH29pUVyqP1BpiUsavg0IID4OII9QBmnHNOz0T109+7dr6H1QCJZU1VkSx0/1Rrr7oGhj2tB0B8ilmoN7NqM/uWmTWa2aCZ1ZnZV8ysZArv8Q4z+5KZPW9m7WbmzOzVy7zGXeLfb6f/kwGYriPN3Tre0iNJyslI0+203gCTVpidocVleZIkJ2n/WdasB/B2Mfn+28yWSNomqULSTyUdlnSNpE9KutvMbnDOtU3irT4u6QFJA5KOS5rsHwSnJT16kccbJvl6ADMoeoLsHasqlJOZ5mE1QOJZW12kk629koIbUV2/pMzjigDEm1g1tX5DwUD/Cefc10YfNLMvS/qUpC9K+ugk3ufvJH1WwT8KFkg6Ncnr1znn/tdUCgYwO5xzY0I9G04BU7dmfpF+tqdRASedae/Thb4hleRmel0WgDgy7fYbM6uVdJekOklfH/f05yT1SvqQmeVd7r2cc6875w4459g2D0gSB5u6wiOMeZlpunUFrTfAVOVlpWtpRX74eF8DLTgAxopFT/3todvnnHOB6Cecc92SXpOUK+naGFxrIsVm9mEz+ysz+7iZzeS1AExB9Cj9navnKjuD1hvgSqyritqI6iyr4AAYKxbtNytCt0cneP6YgiP5yyU9H4PrXcx6Sf8e/YCZ7ZH0Iefcvsm8gZntmOCpldOsDUhZzrkxS1neu5bWG+BKrZ5fqLTdJn/AqbFjQK3dgyoryPK6LABxIhYj9aML6E70XeDo48UTPD9dX5Z0g6RySQWSrpb0hIJB/wUzq5qh6wK4jP1nu3SmvU+SVJCVrpuXl3tcEZC4sjPStGJuQfiY0XoA0WZjnfrR5XXdTLy5c+7TzrltzrlW51yPc267c+5BST+SVCbpM5N8n80X+6fgpF0AV+DpfY3h+++g9QaYtrVjNqKirx5ARCxC/ehvlaIJni8cd95s+Wbo9uZZvi4AvX3Vm3tZ9QaYtlXzCpWRFhwra+ke1LmuAY8rAhAvYhHqj4Rul0/w/LLQ7UQ99zPlfOj2sqvuAIi9PQ2darjQL0kqyE7XTctovQGmKzPdp5XzCsPHextowQEQFItQ/2Lo9i4zG/N+ZlagYL97v6TZ3t11dAWck7N8XQCSntkbab1551XzlJk+G91+QPJbP64Fx7kZ6W4FkGCm/SnrnDsh6TlJixTcETba5xUcKf+uc6539EEzW2lm015Vxsw2XWz9ezNbp+CGV5L0veleB8DUBAK03gAzZdncAmWF/khu7x3S2Y5+jysCEA9itaPsxyRtk/RVM7tD0iFJWyXdpmDbzWfHnX8odGvRD5rZjZI+Ejoc3WVjmZk9OnqOc+4Po17yCUnvNbMXJNVLGlRwCcq7JaVJ+ldJP5zGzwXgCuyq71BjZ7DXtygnQzcuZUt7IFYy0ny6an6hdp4Jtt7sqe9QdUmux1UB8FpMQr1z7oSZbZH0BQUD9T2SmiR9VdLnnXPtk3yrpZIeHvdYxbjH/jDq/pMKTsRdp+AmWNmS2iT9XNK/OueemtpPAiAWokfp775qnjLSaL0BYml9dXE41O9t6NS71lbKZ3aZVwFIZrEaqZdzrl7SI5M896K/eZxzj0p6dArXfFLBYA8gTgQCTs/uo/UGmEm15fnKz0pXz+CIugdHdPJ8r5ZW5F/+hQCSFsNnAGJqx5kL4WX2SnIzdP2SOR5XBCSfNJ+NWbN+dz2r4ACpjlAPIKbGtN6sqVQ6rTfAjNhQHdmo/UBjp4b9AQ+rAeA1Pm0BxIw/4PRMVOvNfbTeADOmuiRHpXmZkqTBkYAOn+v2uCIAXiLUA4iZt+radb57UJJUlp+prYtLPa4ISF5mpg0LIqP1e2jBAVIaoR5AzIxtvZlH6w0ww9ZHteAcae5W/5Dfw2oAeIlPXAAxMeIP6Of7o1a9WTvfw2qA1FBekKWq4hxJwfa3A42dHlcEwCuEegAx8eapdrX2DEkKBo1raL0BZsV6VsEBIEI9gBh5OmqC7D1r5inNx0Y4wGxYV10c3p79VGuvOvuHPa0HgDcI9QCmbcQf0C/2nwsf37uO1htgthTmZGhxeZ4kyUna28BoPZCKCPUApm3biTa19wZbb+YWZmnLwhKPKwJSS/Sa9ayCA6QmQj2AaXtqT2P4/r1r58tH6w0wq9ZUFSk99N9dY+eAWkK7OgNIHYR6ANMyMOzXL6Nab+7fQOsNMNuyM9K0Yl5B+HgPLThAyiHUA5iWl460qHtwRJK0cE7umJU4AMye6DXr9zR0yjnnYTUAZhuhHsC0RLfe3L9+vsxovQG8sGJegbIzgh/r7b1Dqr/Q73FFAGYToR7AFeseGNbzh1rCx/evp/UG8EpGmk9XzY98U8aEWSC1EOoBXLFfHWzW4EhAkrRyXoGWzS24zCsAzKToFpy9ZzvlD9CCA6QKQj2AKzam9YYJsoDnasvzVJCdLknqHRzRifM9HlcEYLYQ6gFckfbeIb16rDV8/G42nAI85zPTuipacIBURKgHcEWe3dekkdBX+5sXlmhBaa7HFQGQpA0LIpu/HWjq0lCoRQ5AciPUA7giT+0eu+oNgPgwvzhbZfmZkqShkYAOn+vyuCIAs4FQD2DKGjv69WZduyTJZ9I9ays9rgjAKDPT+gWRCbO7acEBUgKhHsCUPb03Mkp/w9IylRdkeVgNgPE2RK2Cc7S5W32hDeIAJC9CPYApG7/hFID4Mic/S9UlOZKkgJP2NXZ6XBGAmUaoBzAlx1t6tP9ssEc3M92nd66Z53FFAC4mes16VsEBkh+hHsCU/GRXQ/j+HSsrVJid4WE1ACayrrpIFrpf19an9t4hT+sBMLMI9QAmLRBwenJXpPXmdzZWeVgNgEspyM7Qsrn54eNdZy54WA2AmUaoBzBpb9W162xHvySpODdDt66o8LgiAJeyqSayZv2u+g455zysBsBMItQDmLQnd58N379vXaUy0/kVAsSzVZWFys4I/nfa3jukurY+jysCMFP4RAYwKQPDfj29tyl8TOsNEP8y0nxaWxWZMEsLDpC8CPUAJuWFwy3qHgiudV1Tmjvma30A8WtTTSTU7zvbqWF/wMNqAMwUQj2ASfnJrkjrzXs2VsnMLnE2gHhRU5qrOXmZkqTBkYAONnZ5XBGAmUCoB3BZF3qH9NKRlvAxrTdA4jAzbYward9JCw6QlAj1AC7r6X1NGvYHV83YsKBYi8vyPK4IwFRsXBBplzve0qOu/mEPqwEwEwj1AC7rJzsjG04xSg8knpK8zPAf407SbnaYBZIOoR7AJZ1u69XOM8EAkO4z3beu0uOKAFyJ6MntO89cYM16IMkQ6gFcUvQE2VuWl2tOfpaH1QC4UmvmFyojLTjBvaV7UI0dAx5XBCCWCPUAJuScGxPqf2cTrTdAosrKSNNV84vCxzuYMAskFUI9gAm9eapdp0M7UBZkp+vOVXM9rgjAdES34Oyp72DNeiCJEOoBTOjxHZEJsu9eP1/ZGWkeVgNgumrL81SSmyFJ6h/262ATa9YDyYJQD+CiegdH9Oy+pvDxg5urPawGQCz4zLRpYWS0fkcdLThAsiDUA7ioZ/c1qW/IL0laWpGvDQuKL/MKAIlgc02JRveDPn6+R/XtfZ7WAyA2CPUALiq69ebBzdUys0ucDSBRFOdmatnc/PDx49vrPawGQKwQ6gG8zem2Xr15ql2SlOYzVr0BkszmhaXh+4/vaJA/wJr1QKIj1AN4myeiRulvXV6uioJsD6sBEGurKguUmxmc+N7UOaBXjp33uCIA00WoBzCGP+D0o+jWmy1MkAWSTbrPp41R82QeowUHSHiEegBjvH6iTY2dwZ0mS3IzdPtK1qYHktGWRZEWnF8dbFZbz6CH1QCYLkI9gDEe3xEZsXtgQ5Uy0/k1ASSjuYXZWlCSI0ka9o/dPRpA4uHTGkBYZ/+wfrH/XPiY1hsguW2JmjD72PZ6OceEWSBREeoBhD29t1GDI8Ft41dXFuqq+UUeVwRgJq2tLlJOaKfoo8092l3f4XFFAK4UoR5A2H+9FWm9YZQeSH7ZGWm6d11l+PiHb57xsBoA00GoByBJ2n+2U3sbOiVJmek+/c5G1qYHUsFDW2vC95/a06jOvmEPqwFwpQj1ACRJP4gaobt3baWKczM9rAbAbNm4oFirKgslSQPDAf14V8NlXgEgHsUs1JtZtZl9y8wazWzQzOrM7CtmVjKF93iHmX3JzJ43s3Yzc2b26iRet9rMHjOzFjMbMLMjZvZ5M8uZ3k8FpIbewRH9NGrli9+/puYSZwNIJmamD14b+W/++2+cYcIskIBiEurNbImkHZIekfSmpH+UdFLSJyW9bmZzJvlWH5f0PyRdL2lSa2uZ2VZJb0l6j6RfS/onSV2S/lrSr8wsa/I/CZCafranUb1DfknSkvI8Xb1o0n+LA0gCD2yoUl5oh9njLT1681S7xxUBmKpYjdR/Q1KFpE84597jnPtL59ztCob7FZK+OMn3+TtJayTlS3r35U42szRJ35aUK+n9zrmHnHN/IWmrpB9JukHSp6b6wwCpJnpy3O9fUyMz87AaALMtPytd74maR/O9N5gwCySaaYd6M6uVdJekOklfH/f05yT1SvqQmeVd7r2cc6875w445/yTvPwtklZJ+o1z7qmo9wlI+vPQ4UeNhAJMaP/ZTu2JmiD7vk2segOkog9sXRi+/4v9TWplh1kgocRipP720O1zoTAd5pzrlvSagiPp18bgWhNd+xfjn3DOnZR0VNJCSbUzcG0gKUSP0t+zZp5K8pggC6Si1fMLtammWFJwh9nHtzNhFkgksQj1K0K3Ryd4/ljodnkMrjVj1zazHRf7J2llLAoF4lHv4Ih+ursxfMwEWSC1RY/W/+DN0woEmDALJIpYhPrRLSc7J3h+9PHiGFwrnq4NJLyn9zaqZ3BEklRbnqdrFpde5hUAktm96ypVnJshSapv79dvjp33uCIAkzUb69SP9rN78ef+pK/tnNt8sX+SDs9siYB3fhA1Ge4hJsgCKS87I03vj5pX830mzAIJIxahfnQ0vGiC5wvHnRdLXl4bSGgHGqMmyKb59F4myALQ2B1mnz/UrLMd/R5WA2CyYhHqj4RuJ+pbXxa6najvPVGvDSS0/3j9dPj+3WvmqZQJsgAk1Zbn64alwe1lAm7s7woA8SsWof7F0O1dZjbm/cysQMG14vsl/TYG1xrvhdDt3eOfCC21uVzSaQU3wgIQ0tE3pCd3R/Z3+9B1Cy9xNoBU8/B1i8L3f/jmGfUPTXalaQBemXaod86dkPScpEUK7ggb7fOS8iR91znXO/qgma00s1isKvOypEOSbjaz+6Pe36fgRlaS9E3HftfAGI9tr9fAcHAF2tWVhdqykB1kAUTcsWquFpTmSJI6+4f1k12T2uQdgIfSY/Q+H5O0TdJXzewOBYP2Vkm3Kdj68tlx5x8K3Y6ZlWdmN0r6SOgwP3S7zMweHT3HOfeHUff9ZvaIgiP2T5jZE5LOSLpD0hYF18j/x2n+bEBS8Qec/uO3ka/TH75+IRNkAYyR5jM9fN0i/c0zwY/rR7ed0u9fs4DfFUAci8nqN6HR+i2SHlUwzH9a0hJJX5V0nXOubZJvtVTSw6F/7ws9VhH12MMXufYbkq6W9FMFd7b9lIITZ78g6R3OObbEA6K8eLhF9e3BiW/FuRl6YEPVZV4BIBX97tULlJeZJkk62tyjbScm+1EOwAuxGqmXc65e0iOTPPeif+o75x5V8A+DqV77oKQHp/o6IBV95/W68P3f27JA2RlpntUCIH4VZmfo/Zur9Z3QRNlvv3ZKNywt87gqABOZjXXqAcSJE+d79MqxVkmSz6QPXssEWQAT+4PrF4XvP3+4Rafbeic+GYCnCPVAColemi44ES7Xw2oAxLsl5fm6dUW5JMk56dFtdd4WBGBChHogRfQMjuiJHQ3h4+gl6wBgIo/csDh8//HtDeoeGPawGgATIdQDKeLHOxvUMzgiSVpSnhfeXAYALuXmZWVaUp4n6e2DAwDiB6EeSAHOOX0n6mvzh69fxNJ0ACbFzPSHUaP1j26rkz/A9i9AvCHUAyng5aPndeJ8cIJbfla63rup2uOKACSS922qUmF2cMG80219+uWBcx5XBGA8Qj2QAv7tlVPh+w9uqVZ+VsxWswWQAnIz08eslvXNl0+IzdqB+EKoB5LcgcZOvXo8sozlh6O+RgeAyXrkhsXKTA/Ghr0NnXqdzaiAuEKoB5Jc9Cj9PWsrWcYSwBUpL8jSg5sjrXv//PIJD6sBMB6hHkhijR39+tmexvDxn9xc62E1ABLdn9xcK19ojv0rx1q1/2yntwUBCCPUA0ns0W11GgmtUrF1canWVRd7XBGARLZwTp7etbYyfPwvvznpYTUAohHqgSTVNTCsH7xxJnzMKD2AWPhvtywJ339mb6POtPV5WA2AUYR6IEn915v1Yzabum1FhccVAUgGa6qKdOPSMklSwEn/+gqj9UA8INQDSWjYH9C3XotMkP3jm2rl87HZFIDY+GjUaP1j2+vV2jPoYTUAJEI9kJSe3dekps4BSVJZfqbes7HK44oAJJMbls7RmqpCSdLgSGDMjtUAvEGoB5KMc07/8nLk6/CHr1uk7Iw0DysCkGzMbMxo/aOv1amzb9jDigAQ6oEk89KR8zrY1CVJys7wjdkFEgBi5V1rKlVbnidJ6h4c0b9HtfwBmH2EeiCJOOf01ReOhY8fumahSvIyPawIQLJK85k+ecey8PG3Xz3FaD3gIUI9kEReP9GmXWc6JEmZaT6WsQQwo+5bN19Lokbr/+1VVsIBvEKoB5LI1144Hr7/4JZqzSvK9rAaAMkuzWf6RPRo/Wt16ugb8rAiIHUR6oEksb2uXa+fbJMU/KCNnsQGADPlvnXztbQiX5LUMzjCuvWARwj1QJKIHqX/nY1VWlCa62E1AFLF+N76R1+rU3svo/XAbCPUA0lgb0OHXj56XpLkM+ljtzJKD2D23Lu2UsvnBkfre4f8+jdG64FZR6gHksD/jhqlv2/dfNWW53tYDYBU4/OZPnnH8vDxd7YxWg/MNkI9kOAONXXpuYPN4eOP37bUw2oApKp3rZmnFXMLJAVH6//5peOXeQWAWCLUAwkuepT+nVfN1Yp5BR5WAyBV+Xym/35npLf+O9tOq769z8OKgNRCqAcS2IHGTj2zryl8/Ge3LbvE2QAws+5eM08bFhRLkob8AX3puSMeVwSkDkI9kMC+/NzR8P27Vs/V2uoiD6sBkOrMTH91z6rw8ZO7G7X/bKeHFQGpg1APJKgdpy/o+cMtkiQz6dN3rfC4IgCQrllcqnesnhs+/n+ePSTnnIcVAamBUA8kqH/4ZeRr7fvXz6eXHkDc+Iu7VyrNZ5KkbSfa9FJoyV0AM4dQDySg1463jtk99lN3Lr/MKwBg9iytyNf/cfWC8PHfPntY/gCj9cBMItQDCcY5p7+PGqX/3S3VWlSW52FFAPB2n7xzmXIz0yRJR5q79aMdDR5XBCQ3Qj2QYH59qEW76zskSZlpPv2ft7PiDYD4U1GQrT+9ObK79Zd+dUR9QyMeVgQkN0I9kEACATdmibgPXFuj+cU5HlYEABP7yE2LVV6QJUlq7hrUP790wuOKgORFqAcSyM/2NurwuW5JUm5mmj52K7vHAohfeVnp+sxdkTk///LySZ1q7fWwIiB5EeqBBDEw7B/TS//IDYvCI2AAEK8e3LxA66M2pPrcUwdY4hKYAYR6IEE8uq1ODRf6JUkluRn6k6heVQCIVz6f6W8eWCMLrnCp3xw9r18eOOdtUUASItQDCaCtZ1Bff+F4+Pi/37lcRTkZHlYEAJO3trpIH9y6MHz8hZ8dZNIsEGOEeiAB/NPzx9Q9GPwArC3P00NbazyuCACm5jN3rdCcvExJUmPngL4WNVABYPoI9UCcO97So++/cSZ8/FfvWqWMNP7TBZBYinIz9JfvWhk+/rdXTup4S4+HFQHJhWQAxLn/99lD4Z0Yr6udoztWVXhcEQBcmfdtqtaWhSWSpGG/01//dD+TZoEYIdQDcWzb8VY9f7hFkmQmffbeVbLR2WYAkGB8PtMXHlgjX+jX2LYTbXqcnWaBmCDUA3HKH3D6m2cOhY/ft6laa6qKPKwIAKZv9fxCPXLD4vDx//2zg2rq7PewIiA5pHtdAICLe3x7vQ42dUmSMtJMS8rz9YOo3noASFSfuWuFnj/UrLq2PnUPjugvf7RPjz5yNd9EAtPASD0Qhy70DunvfnE4fHzTsnKWsASQNHIy0/T3D64Pr13/8tHzenw7bTjAdBDqgTj0988d0YW+YUnBjaZuXlbucUUAEFtXLyrVI9dHteE8fVCNHbThAFeKUA/EmT31Hfrhm5E2m/vWzVdmOv+pAkg+//OdK7S4LE+Sgm04P97HajjAFSIpAHHEH3D6v366X6OfaSvmFmjlvAJviwKAGZKTmaa/f/+6cBvOb46e12Pb670tCkhQhHogjvznW2e0t6FTkpSZ/v+3d9/xVdX3H8df35s9SICEbUwe1wAAIABJREFUhBH2CrKE4ABUBC0O3MVZqVotilgctdpfh1Xb2mEVZ6vWgYgWqzjqXgxliQxBZI8wQ0ISCNnz+/vj3FxCTFi5uSc3eT8fj/s4nO8553s+N/fk8sn3fM/36+GCQR304JiINGvDurXlZzVGw3nwvTVs2atJqUSOlZJ6kSYit7CMv3+83rc+aVRPEmIjXIxIRCQw7h7blx7ebjiFZZXc9toKSsorXY5KJLgoqRdpIv7+8Tryip2HYzu3jWLSmT1djkhEJDCiwkN44uohhIc4acmajAM89OHaIxwlIjUpqRdpAr5Jz2XmNwf7kd5/YX8iw0JcjEhEJLAGdIrndxf0861PX7SNj77LcDEikeCipF7EZSXlldz75irf+tn9kjirX7KLEYmIuGPCqV05t3973/o9s1axI7fIxYhEgoffknpjTIox5kVjzG5jTKkxJt0Y85gxps0x1tPWe1y6t57d3npT6tk/3Rhj63nt8c+7E2k8T87eyJbsQgBiI0J58OIBLkckIuIOYwx/Gz+IlDZRAOSXVHDbf1ZQVlHlcmQiTV+oPyoxxvQEFgJJwLvAOuBk4HbgXGPMSGttzlHUk+Ctpw8wG5gJpAI3AOOMMcOttVvqODQPeKyOcj0+L03a97vzeGbewUv63vNS6dg6ysWIRETcFR8VxlPXDGX8vxZSUWVZuWM/D324lvsv6u92aCJNml+SeuCfOAn9FGvtk9WFxphHgTuBPwO3HEU9D+Ek9FOttXfVqGcK8Lj3POfWcdx+a+39xx29iAsqKqu4d9YqKqucQelP7taWn5zcxeWoRETcd2Ln1vz6vFT+9IHzsOy0hemktm/FVfqOFKlXg7vfGGN6AGOBdODpWpv/ABQCE4wxMUeoJwaY4N3/D7U2P+Wt/xzv+USC3vPzt7J61wHAGZP+rz8eiMejMelFRABuPK37If3rf//uar7ecsSb/iItlj/61I/xLj+11h7S6c1amw8sAKKBU49Qz3AgCljgPa5mPVXAp97V0XUcG2GMudYY8xtjzO3GmNHGmGMaOsQYs6yuF073HxG/2ppdyNTPNvjW7zi7Nz3axboYkYhI02KM4dErB9OvQxwA5ZWWSa8u14OzIvXwR1Lf17vcUM/2jd5ln0aspz3wCk43n8dw+uNvNMaMOsI5RQKuqspy76xVlHof/DqhQxw/P103oEREaosOD+X564aRGBsOOJP03fTyUgpKK1yOTKTp8UdSH+9d5tWzvbq8dSPV8xJwFk5iHwMMBJ4FugEfGWMGH+G8AFhr0+p64Tz0K+I3L8zfypKtuQCEeAx/Hz+IsBCNLisiUpdOraN4dkKab2Kq9Zn53DFzhe95JBFxBCKTqO4k3NDfvjrrsdY+YK2dba3NtNYWWWtXW2tvAR7F6c5zfwPPK+I36/Yc4OFP1vvWbz2zJwM6xR/mCBERSevalocuG+hb/3xtFn98fw3WKrEXqeaPpL66Bb2+zCSu1n6NXU+1Z7zLM45yf5FGVVpRyR0zv6Ws0ul2M7BTPFPO6u1yVCIiwWF8Wgo3n3Gwq+K0hek8/sXGwxwh0rL4I6mvbnasr898ddZSX195f9dTLcu7POyoOyKB8uhnG1i3x3kGPCLUw9QrT1S3GxGRY3DPuamMG9TBt/7Y5xuZtmCrixGJNB3+yCjmeJdjjTGH1GeMaQWMBIqBxUeoZ7F3v5He42rW48EZNrPm+Y5kuHdZ12RVIgH19ZYcnvvy4KX4m/P70StJo92IiByLEI9h6hUncnrvRF/Z/e+t4Z0Vu1yMSqRpaHBSb63djDPcZDdgcq3ND+C0lE+31hZWFxpjUo0xhwwVaa0twBnBJoYf9oO/zVv/JzVnlDXG9DfGtK0dkzGmK87Y9gAzjvlNifhRfkk5d/13JdVdP0/vnciEU7u6G5SISJAKD/Xw7IQ0hnQ5OG7G3W+sZPa6TBejEnGfv+7934rT3eUJY8w7xpi/GGNm48wmuwH4ba3913pftf3Gu/9dxpgvvPW8gzObbBY//KPhcmC3MeYjY8w/jTF/M8a8iTNiTS/gQ+AffnqPIsfMWst9737Prv3FgDP9+cPjB2uSKRGRBogOD+Wl60+ib7JzY7+iyjJpxnK+2rjX5chE3OOXpN7bWj8MmAacAvwS6Ak8AQy31h7VFHDe/YZ7j+vlrecUnGEr07znqWkO8DbQHbgGuAsYBcwHrgMusNaWNeS9iTTEG8t28naN28J/vnQA7eMjXYxIRKR5aB0dzvQbTyalTRQApRVV3PjyUuasyzrCkSLNU6i/KrLW7gBuOMp9622mtNbmArd7X0eqZx4w72hjFAmkDZn53Pfuat/6+LQULhjU0cWIRESal+S4SF696RSufm4xu/NKKKuoYuIrS3nqmqGc07+92+GJBJSG3hBpBEVlFdz66nJKyp3hK3snxfLgxf1djkpEpPnpmhDD6zcPp3Nbp8W+vNJy66vLeX/VbpcjEwksJfUijeC+d79nU1YBAJFhHp7+yVCiw/12Y0xERGro3Daa1ycOp3uiM4p1ZZVlyn9W8NbynS5HJhI4SupF/GzWsp28uezgfyQPXjyAPsmtDnOEiIg0VMfWUbw+8VTfcMFVFu7670qembdZM89Ki6CkXsSPNmbm87t3Dvajv2xIJy5PS3ExIhGRliMpLpKZE08ltf3BhpS/frSO37+7mgrvbN4izZWSehE/ySsuZ+IryygurwSgZ7sY/njJAIzR8JUiIoGSGBvB6xOHc0r3g9PYzFi8nZtfWUZRWYWLkYk0LiX1In5QWWW5Y+YKtmY7c6xFhYXw9E+GEhOhfvQiIoEWHx3G9BtP5qLBB0cc+2JdFlc9t5is/BIXIxNpPErqRfxg6mcbmLP+4KQnD18+iNT2cS5GJCLSskWEhvDYlScy6cyevrJVO/O4+KkFrNq538XIRBqHknqRBvrwuwyemrPJtz7pzJ4aj15EpAnweAz3npvKny4ZQPVE3hl5JYx/ZhFvLN3hbnAifqakXqQB1u05wN1vrPStj+rTjrvH9nUxIhERqe3aU7vywvUnERfpdIksq6jiV2+u4r53V1NWoQdopXlQUi9ynHILy5g4fRlFZc6Dsd0SonniqiGEePRgrIhIUzO6bxL/u+00+iTH+sqmL9rGT55fTNYB9bOX4KekXuQ4lJRXMnH6UrbnFgEQEx7Ccz8dRnx0mMuRiYhIfbolxvD2rSM5f2B7X9k36fs47/GvmLs+y8XIRBpOSb3IMaqqstz9xkqWbtsHgDHw6JUnaoIpEZEgEBMRytPXDOXX56X6+tnnFJZx/Uvf8NCHa9UdR4KWknqRY/Twp+t5f1WGb/235/fjnP7tD3OEiIg0JcYYbhnVkxk3nkK7VhG+8ue+3MLlzyxkW06hi9GJHB8l9SLH4LWvt/OvuZt969cN78qNp3V3MSIRETleI3ol8tHtp3Nm33a+spU78xj3xHz+u3QH1loXoxM5NkrqRY7S3PVZ/P7d1b71s/slcd+F/TVjrIhIEEuMjeDF607id+P6ERbifJ8XlFZwz5ur+Pn0pZqsSoKGknqRo/Dtjv1MfnU5lVVOq83ATvE8cbVGuhERaQ48HsNNp/dg1qQRdEuI9pV/vjaLsVO/5P1Vu12MTuToaA57kSNYvyef615cQqF36MpOraN44bphRIfr10dEgt9rX293O4Qm5foR3fn4+z0s3pIDwP6icm57bQUfr97DAxf1JyE24gg1iLhDLfUih7Etp5BrX/iavOJyANrGhPPyz04iKS7S5chERKQxhId6uGhwR1696RQ6xh/8rn9/VQZnPzqPt1fsVF97aZKU1IvUY09eCT95/mv25pcCEBsRyss3nEyvJA1dKSLS3I3slcjHd57B5WkpvrJ9ReXc+fpKrnvpG3Z45ykRaSqU1IvUIaeglGtf+Jqd+4oBiAj18MJ1wxiYEu9yZCIiEihxkWE8fPlgpt1wEp1aR/nKv9ywl7FTv+T5r7ZQUalx7aVpUFIvUktuYRkTXljCpqwCAEI9hmeuTeOUHgkuRyYiIm44s28Sn955BjeM7Eb1gGfF5ZX86YO1jHtiPl97+9+LuElJvUgN2QWlXPPvxazJOAA4s8VOvfJERqcmuRyZiIi4KSYilD9c2J+3Jo2gb40ZxNdn5nPlc4u58/VvyTqg4S/FPUrqRbyy8ku4+rnFrNuTDzgJ/d9+PIgLB3d0OTIREWkqhnRpw3u/OI3/Oy+V6PAQX/nbK3Yx5pF5vDB/q7rkiCuU1IsAmQdKuOq5xWz0drnxGHjk8sFcMayzy5GJiEhTEx7q4eZRPfnil6MYN6iDr7ygtII/vr9GXXLEFUrqpcXLyCvmqucWs2VvIQAhHsPUK0/ksqEpRzhSRERasg7xUTx9zVBevekUeraL8ZWrS464QUm9tGgbM/O57J8L2ZrtJPShHsOTVw/h4hM7uRyZiIgEi5G9Evno9jP4dT1dcp6Zt5mS8koXI5SWQEm9tFjfpOcy/plFZOQ5rShhIYanfzKU8wd2OMKRIiIihwoP9XBLPV1y/vrROs5+dB7vrdytiauk0Siplxbpk+/3cO3zB2eKjQkP4YXrTuKc/u1djkxERIJZdZecGTce2iVn575ifvGfFVz2r4Us27bPxQiluVJSLy3OjMXbmDRjGaUVzugEibHhvH7zcM7o087lyEREpLk4rXciH99xBg9e3J820WG+8hXb9/Pjfy1k8mvLNSut+JWSemkxKiqr+PMHa/jdO6up8t797JYQzaxJIxjQSTPFioiIf4WFePjp8G7M/dVobj6jB+EhB9OuD1ZlcNYj83jow7W+u8YiDaGkXlqEvKJybpj2Df/+aquvbFBKPG9OGkHXhJjDHCkiItIw8VFh/N/5/X7Q376ssornvtzCmQ/PYdqCrZRW6GFaOX5K6qXZ25SVzyX/XMBXG7N9ZWf3S+Y/Pz+VxNgIFyMTEZGWpHPbaJ6+ZiizJg1nSJfWvvJ9ReXc/94axvxjHm8s3aHJq+S4KKmXZm32ukwuffrgkJUAvxjTi+cmpBETEepiZCIi0lKldW3LW5NG8OTVQ0hpE+Ur37W/mF+9uYpzHvuSD7/LoKpKI+XI0VNWI81SeWUVj3y6gWfmbfaVRYWF8I/LBx9y61NERMQNxhguHNyRH52QzIzF2/jn3M3kFpYBsHlvIbe+upwBneK4e2xfRvVphzHG5YilqVNLvTQ7O3KLuOLZRYck9J1aR/HmpOFK6EVEpEmJDAvhptN78OU9o7nz7D7E1riLvHrXAa5/6RuufHYxS7bmuhilBAMl9dKsfPhdBuc/8RUrtu/3lZ3Rpx3/u20k/TtqhBsREWmaYiNCuf3s3nx1jzNSTkTowRRtSXouVzy7iOtfWsK3O/YfphZpyYxmNjs8Y8yyoUOHDl22bJnbochh5JeU89CHa/nPkh2+slCP4Vfn9OXnp/fA4wm+25bdfv0BAA9dOrBRz/Obt78LyHkaw7HGXr1/zWPqq6Ouff0Rw9HW4dbnUvN9H+/5a8ceiPdyuLhrb/NHLA39zAL5+Tal6ytYVP980v86zrUYMg+U8OTsjcxcsoOKWn3rT+uVyOTRvTi1R1t1y2lm0tLSWL58+XJrbdqxHqs+9RL05qzP4jdvfUdGXomvLKVNFE9ePYQhXdq4GJmIiMjxSY6L5E+XDGTi6T157PMNvP3tLqrbYedvymb+pmzSurZh8uiejO6bpORelNRL8NpfVMYf31/LrOU7DykfN7ADD102kPiosHqOFBERCQ5dEqJ59MoTueXMnjw9ZxPvrdztm0Bx2bZ9/GzaUvp1iGPy6J6cN6ADIUF4Z1r8Q0m9BB1rLe+vyuDB99ewN7/UV942JpwHLurPBYM6qMVCRESalT7JrXj8qiHc9aM+PDNvM28u20l5pZPdr804wG2vraBHuw1MGtWTS4Z0IixEj022NErqJah8vzuPB95b84NRAC4a3JE/XHgCCZpMSkREmrGuCTH85bJBTDmrN//+ciuvLdlGSbkzWdWWvYX86s1VTP1sA9eN6MZVJ3fRXesWREm9BIWcglIe+WwDM5dsp+bzQkmtIvjTJQMY27+9e8GJiIgEWIf4KO678AQmj+7JSwvSeXlROvklFQDszivhLx+t4/EvNjI+LYXrR3SjR7tYdwOWRqekXpq0gtIKXl6YzrPzNnPA+2UFzsg214/oxi/O6q1WCBERabESYiO4+5y+TBzVg1cWbePF+VvJ8U5iVVRWyfRF25i+aBtn9m3Htad0ZXRqkvrdN1NK6qVJKiqrYPqibTw7bzP7isoP2TaqTzt+f8EJ9EpSq4OIiAhAXGQYk0f34sbTuvO/b3fz4oKtrNuT79s+d/1e5q7fS8f4SK46uQtXntSZ5LhIFyMWf1NSL01Kfkk5M5fs4NkvN5NdUHbItu6JMfxuXD/GpGroLhERkbpEhoVwxUmduXxYCos25/Digq18sS7LNxzm7rwSHv1sA49/sZHRfZMYn9aJManJhIfqwdpgp6RemoQduUW8tCCd/y7dQUFpxSHbUtpEMWVMby4dqqf5RUREjoYxhhG9EhnRK5HtOUW8tmQ7byzd4euaU1ll+XxtJp+vzaRNdBgXDe7Ij9NSGNgpXg1nQUpJvbimqsry9dZcpi3cymdrMqk1YR4d4yO5bUxvxqelqAVBRETkOHVJiObX56Vy549688n3mcxYvO2QUeT2FZXz8qJtvLxoG90TYxg3sAPjBnUgtX0rJfhBREm9BNyO3CJmLd/JrOU72ZFb/IPtPdvF8LPTujM+LYWI0BAXIhQREWl+IkJDuGhwRy4a3JGt2YW8tXwns5btZHeNGdm3Zhfy1JxNPDVnEz3aOQn+j05IZkDHeDx6wLZJU1IvAZFdUMrnazJ599vdLNqSU+c+p/dO5MbTunNG73b64hAREWlE3RNj+OXYvtx5dh8Wb8nhzeU7+WT1HgrLKn37bNlbyJOzN/Hk7E0ktYpgTGoSZ/VL5rReiUSFq9GtqVFSL41m1/5iPlm9h4+/38PS9NwfdK8BiIsM5eITOzFheFf6JLcKfJAiIiItmMdzsO99yaWVzNuwl/dXZfDF2kyKaiT4WfmlzPxmBzO/2UF4iIe0rm0Y2SuBkb0SGdgpnlA98+Y6JfXiN/kl5Xy9JZf5m7JZsCmbjVkFde7nMXB673ZcPiyFs/slExmmv/ZFRETcFhkWwjn923NO//YUl1Uyd30Wn63JZM76rEOGly6rrGLRlhwWbcnhH59uoFVEKCd1b0ta1zYM69qGwZ1b6/92Fyipl+OWkVfM8m37WbF9H8u272PVzjwq62qOB4yBYV3bcE7/9lwwqCPt4zU2roiISFMVFR7CeQM7cN7ADlRWWVZs38fna7OYvS6TDZmHNtrll1Ywe10Ws9dlAc4Ekf07xXNiSjz9O8UzoGM8vZNjNYJdI/NbUm+MSQEeBM4FEoAM4B3gAWvtvmOopy1wH3AJ0AHIAT4G7rPW7mzMc0vdKqss23OLWJdxgLV78lmXcYBVO/PYc6DksMeFh3g4pUdbzhvQgbNPSCKplRJ5ERGRYBPiMQzr1pZh3dry6/NSyTxQwsLN2czfmMOCTdk/yAcqqiwrd+xn5Y79vrLwEA9927eiT3IreifH0ic5lt5JrejUOkrP0fmJX5J6Y0xPYCGQBLwLrANOBm4HzjXGjLTW1v105KH1JHjr6QPMBmYCqcANwDhjzHBr7ZbGOHdLV1FZRVZ+Kbv2F7M1u5BtOYWkZxexNbuQrdmFFJdXHrkSoH/HOE7rlcjIXomc1K2tHqQRERFpZpLjIrl0SAqXDknBWsvW7EKWbtvHsvR9LN2Wy+a9hT84pqyyiu925fHdrrxDyiNCPXRuG03XttF0SXCWKW2iSYqLIDkukoSYcPXXP0r+aqn/J05SPcVa+2R1oTHmUeBO4M/ALUdRz0M4Cf1Ua+1dNeqZAjzuPc+5jXTuZqWisoqC0grySyp8y9zCMnIKS8kpKCO7oJTsglIy8krI2F9CVn5JnQ+yHk50eAiDU1ozpEtrhnRpw9AurUmIjWicNyQiIiJNjjGGHu1i6dEuliuGdQZgX2EZ3+7cz/e78li96wDf7cpj1/4fDmENUFpRxaasAjYd5jm8hNgIkuMiSG4VSVJcJEmtIkiKiyA+Koy4yDBaRYYSF+VdRoa12P78DU7qjTE9gLFAOvB0rc1/ACYCE4wxv7TW/vBPt4P1xAATgELvcTU9hZOgn2OM6VHdWu+vczdFr3+znT15pZRXVlFeWUWZd1leYQ9dr7SUlFdSUFpBQUkF+d7l0basH612rSJIbd/K+4qjX4c4+iTH6q9nEREROUSbmHBG901idN8kX9m+wjLW7jnApqwCNmYWsCEzn41ZBeR6Z7itT5WFvfml7M0vZTUHjur84aEe4rwJfnRECBGhIYSHeIgI83iXtdZDPYR4DCEegzEGj4EQY4gI8zDxjJ4N+lkEkj9a6sd4l59aa6tqbrDW5htjFuAk3qcCXxymnuFAlLee/Fr1VBljPsVJ0kcD1V1w/HXuJmf6om18v/voLl5/SYyNoGPrSLq0jaZ7YgzdEmLolhhN98RY2saEBzQWERERaT7axIQzomciI3omHlKeV1zO9pwituUWsi2niO05Rew5UELmgRL25peSc4Skvy5lFVVkF5SRXXDsx9bUKiI0qJJ6Y+0x9rmoXYExDwN3A3dbax+pY/tTwGTgVmvtvw5Tz2ScFvmnrLW/qGP73cDDwN+ttff689zefZfVs2lwVFRUSL9+/Q53uN9t3ltAcVnDWts9HkOIMb5liMcQ6jGEhhhCPR5CQgxhHkNYiIewEA+aCbppWe3td9ixdVSjnme395ZoY5+nMRxr7Ltr3P6tPqa+Oura1x8xHG0dbn0uu2vdIj+e89eOPRDv5XBx197mj1ga+pkF8vNtStdXsKj++QzoFO9yJGItVFQ5PRMqKqsor3KWFZWWSmuprLJUeZfV/25gauvj8RhO6BDnn8qO0tq1aykuLs611iYc67H+aKmvvuLz6tleXd66Eerx17kPp7K4uDhv+fLl6Q2oI5BSvct1rkYhfpOeeVyHHfN1cJznaRKOJ/baxxyujqOp3x8/v7rqaGC9Df4+aMj5j+Vn7G9HOpe/YmnoZxagn0mq91w/uA6C+fc+EJY3r5+P8oPjsDwj4KfsBkfZz6iWQIxTX93+29C/m46nnqM+xlqbdswRNUHVdxyay/uR46PrQEDXgTh0HQjoOmgJ/PGUY3VreH33qOJq7efPevx1bhERERGRoOWPpH69d9mnnu29vcsNjVCPv84tIiIiIhK0/JHUz/EuxxpjDqnPGNMKGAkUA4uPUM9i734jvcfVrMeDM4pNzfP589wiIiIiIkGrwUm9tXYz8ClOx/7JtTY/AMQA02uOE2+MSTXGpNbc0VpbALzi3f/+WvXc5q3/k5ozyh7PuUVEREREmpsGD2kJYIzpCSzEmdn1XWAtcArOmPIbgBHW2pwa+1sAa62pVU+Ct54+wGxgCdAPuBjI8tazuSHnbu70IIyArgNx6DoQ0HUgDl0HzZ9fknoAY0xn4EHgXCAByADeAR6w1ubW2rfOpN67rS3ObLCXAB2AHOAj4D5r7c6GnltEREREpLnxW1IvIiIiIiLu8MeDsiIiIiIi4iIl9SIiIiIiQU5JvYiIiIhIkFNSLyIiIiIS5JTUi4iIiIgEOSX1IiIiIiJBTkl9C2OMecEYY72vXm7HI43PGNPbGHOvMWa2MWaHMabMGJNpjHnXGDPa7fjE/4wxKcaYF40xu40xpcaYdGPMY8aYNm7HJo3PGJNgjLnJGPO2MWaTMabYGJNnjJlvjLnRGKP/+1swY8yEGnnATW7HI/6jcepbEGPMhcD/gAIgFuhtrd3kblTS2IwxM4ErgTXAfCAX6AtcBIQAt1trn3AvQvGnOmbZXgecjDPL9npgZEuaZbslMsbcAvwLZyLGOcB2IBm4DIgHZgGXWyUALY53ss7vcL77Y4GfW2ufdzcq8Rcl9S2EMaYdzi/yXKA9MAol9S2CMeZ6YKW1dkWt8lHAZ4AFullrM1wIT/zMGPMJMBaYYq19skb5o8CdwLPW2lvcik8anzFmDBADfGCtrapR3h5YAnQGxltrZ7kUorjAGGNwvvO7A28Bd6OkvlnRLbiW4znvcrKrUUjAWWun1U7oveXzcP7ICwdGBDou8T9jTA+chD4deLrW5j8AhcAEY0xMgEOTALLWzrbWvlczofeW7wGe8a6eGfDAxG1TgDHADTjfBdLMKKlvAbwttZcAt+i2u9RS7l1WuBqF+MsY7/LTOhK6fGABEA2cGujApMnQ73wLZIzpB/wVeNxa+6Xb8UjjUFLfzBljugKPAzOste+4HY80Hd5r4yygCNCXfPPQ17vcUM/2jd5lnwDEIk2MMSYU+Kl39WM3Y5HA8X7ur+A8W/Ebl8ORRhTqdgDSeLwjHLyM82DsFJfDkSbEGBMBvApEAPdYa/e5HJL4R7x3mVfP9ury1gGIRZqevwIDgA+ttZ+4HYwEzH3AEOA0a22x28FI41FLfRPnHYrOHsNrRo3D78R5IPbnStqCWwOvg9p1heC02owEXgf+Eaj3Ia4z3qVGSGhhjDFTgF/ijIY0weVwJECMMSfjtM4/Yq1d5HY80rjUUt/0bQZKjmH/3eCMTQ78GXjJWvthYwQmAXVc10Ft3oR+BnA58F/gWg1r16xUt8TH17M9rtZ+0gIYYybjdMNcA5xlrc11OSQJgBrdbjYAv3c5HAkAJfVNnLX2rOM8tD9O14objDE31LPPRmeEKy5Vf/umrQHXgY/3C/41nIT+NeCn1trKhtYrTcp677K+PvO9vcv6+txLM2OMuQOYCqzGSeizXA5JAieWg98FJd7/72v7tzHm3zgP0N4RsMikUSipb77SgRfq2TYOZ6z6N4AD3n2lGTPGhOO0zF+6EDUDAAABoElEQVQMTAduqD06ijQLc7zLscYYT60xylvhdLkqBha7EZwEljHmXpx+9N8CP7LWZrsckgRWKfXnAUNx+tnPx2kMUNecZkCTT7VAxpi5aPKpFsP7UOxbwPk4X/ATldA3X5p8SgCMMb8HHgSWAWPV5UZqMsbcjzN3hSafakbUUi/S/D2Dk9BnA7uA++q4DTvXWjs3wHFJ47gVWAg8YYw5C1gLnAKMxul281sXY5MAMMZch5PQVwJfAVPq+J1Pt9ZOC3BoItKIlNSLNH/dvctEnKHN6jO38UORxmat3WyMGYaT1J2L8wddBvAE8IBabFuE6t/5EKC+ftLzgGkBiUZEAkLdb0REREREgpzGqRcRERERCXJK6kVEREREgpySehERERGRIKekXkREREQkyCmpFxEREREJckrqRURERESCnJJ6EREREZEgp6ReRERERCTIKakXEREREQlySupFRERERIKcknoRERERkSCnpF5EREREJMgpqRcRERERCXJK6kVEREREgpySehERERGRIKekXkREREQkyCmpFxEREREJcv8Phk7gkGiZT2IAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 248,
       "width": 378
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "seed(314)\n",
    "N = 30\n",
    "mu, tau = 0, 0.7\n",
    "sigma = np.sqrt(1 / tau)\n",
    "D = normal(loc=mu, scale=sigma, size=N)\n",
    "\n",
    "sns.distplot(D, rug=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The true posterior is given by a Normal-Gamma Distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuMAAAI0CAYAAAC3RiczAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAWJQAAFiUBSVIk8AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzde9gkZX3n//cXEA+cRdFNUJCTw8ZNFIyIGCOgSEyixuhuftklwkYTfxLBU06aCJg10c2lICox0dVRk/3FBJOYbBLFIPHAQbIYY1RGdHAUggQE5SQMCN/fH1XN9PT06Xn6UHdVv1/X1VfN01V9V1X3M1Wf/j533RWZiSRJkqTl26npDZAkSZJWlWFckiRJaohhXJIkSWqIYVySJElqiGFckiRJaohhXJIkSWqIYVySJElqiGFckiRJaohhXJIkSWqIYVySJElqiGFckiRJaohhXJIkSWqIYVySJElqiGFckkaIiC0RkRGxseltkSR1k2FcaomIeHodDHuPj0zxmjP7lt+wjO3UZBHxgxHxyoj4PxGxOSJuiYi7I+LbEfEvEfH+iHhRROzV9LZKkhZrl6Y3QNK6PScinpSZlze9IZpORDwEeBPwS8ADhyyyb/34YeAXgK0R8W7gDZl549I2VJK0NIZxqd3+B3BC0xuhySJif+D/AD9SP3Ub8CHgH4Fr6p/3BR4DHA88G9gD+BXgs8AfL3eLJUnLYBiX2unbwMOAZ0bE0zLzU01vkEaLiAcBf822IP6/gZdn5s0jXvLuiNgbeDnwm0vYRElSQ+wzLrXTO4A76n//jyY3RFP5HeAJ9b/fn5n/dUwQByAzv5uZvwMcAVy16A2UJDXDMC61079TBXKAH4uIdXdViYiHRMRrIuIz9QWEWyPiW/XFhT8fETHmtRvri0O31D8/MiJ+NyK+WF+UmBHxvHregX0Xk55cP/f8iLggIm6IiDvqixdfHhEP6FtH1Nvxj/Vy34uIz0XESyds224R8V8i4j0R8fl6e+6JiBsj4pP1Pu++3vdtWhHxUOD/rX/8d6puJ1PLzE2jrguIiMdExKsi4q8j4usRcWf9+EZE/FlEPHvCtv23vs/kkIh4QES8IiL+KSK+GxHfiYhPRcRPD7xu9/r9++eIuLV+by+MiGeOWddT+9b1jPpzPSUiPh0RN9Xt/FN94Wr0ve4BEfHLEXFZRNwcEbdHxKUR8XMT9m3viDg5Ij4YEV+KiNvqi2T/vf6d++WIGNZvX5KWKzN9+PDRggfwdCDrx0uBhwK31D9/dsRrzux7zYYh8/8T8M2+ZYY9PgXsM6L9jfUyW4CjgBuGvP559bIH9j13MnDemHV+GNiZ6iLHPx+z3B+Neb/+ccJ+JXD1sPelr40t9XIbZ/jcXtq3vv8xx9+HR02xfwn8CfCAEW38t77lngBcPKadX6tfsz/whRHL3Ae8aMS6ntq33InAX45Z1zvr1+wFfGLMcq8f8/58bYr35vPADzT9f9uHDx+r/bAyLrVUVt0czq5/fFJEPGctr4+IHwAuogp1UPVjfjbwRODnqIIZwI8BfxsRO49pbneqAP0Q4PeAY4EfpRoRZMuQ5V9KVS3+O+D5wJHA86guVKR+7hTg94EX1Nv2U/VyPwdsqpd7SUScOGKbdgH+FXgj8DNUXxaeDPwX4E+pguNjgL+q+3QvytP7/v33c2x3Z+Ae4G+B04FnUHVpOZ7qvf1ivdzPU30pm+QPqN6jdwHPpHqv/zvwb/X834uIH6b6nA8B/ifV5/xE4DTgu0AA74iIR05Y15lUn/f/x7bP9b8AV9bzX1Z/ru8Fngb8IdWFykdS/V70tumMepuG2Rm4HHg98NNUv4/HAP+V6vcOqj78Hxr3FxZJWrimvw348OFjugcDlfH6uT2Bm9hW5YuB15zZ95oNA/M+1DfvtCHr24kqtPaWefmQZTb2zb8DOGLM9h/I9lXJs4cs8xC2VaNvpArMpw9Z7pHArfVyHxmxvkMnvJ/PAO6t2/jFEcv0tmXjDJ/bl/v2+cFz/H14CGOqulTB+N19n80Of91g+8p4Ai8YsswPAd+v598AbAWeMmS5E/ra+dUh8586sK5XD1lmP6pQ31vXNNv0znV+/if1bcuz5vW5+PDhw8daH1bGpRbLzFupqsdQVfleOM3rIuI/UFWfAT6VmecOafs+4JeB3oWGp01o9n9m5uemWT/VUH6/NmSd3wPeX//4MKruN28bstz1VN0coKrc7yAzvzpuAzLzH6hGOIGqSrsoD6unt2fmnaMWioh9IuJxIx6PHVw+M7+XmdeNai8zE/hVqur5Q5g8BOaHM/P8Ie18Cfhk/ePDgXMz85Ihy10AbK5/HPqZ9Pm/mfmWIW3cQFV5763r/Cm2ab2f/weB/1v/uMjPX5LGcmhDqf3OBV4BPAI4KyI+nJn3TnjNsWz7//+eUQtl5i0R8WdU3UoOiYgDM3PLiMXXMg72X2TmPSPm/Uvfvz80po3ecvtExN6Z+d1xK4yIhwN7s/3Ndno30vmRHV8xN3vU09snLPezVJXsYf6Nqq/2SBGxC9VfDPag6qLR823gP1B3yRjTxJ+OmfcvwHH1v8e18QXgYOCgcds6oY21fP7HTbEuACLiEVR90Hfte7r3ZWaRn78kjWUYl1ouM78XEb8HnANsoOp68P7xr+Jxff++bMKyl1GFcagu+NwyZJnbM3PzkOdHGTdUX3+onna5PQZ+BiAijqGq6D+D6oLXUR42Zt4OohoDfGQ4zswv9v14G/Agqn71c1UH8JdQ9c0/gu2D5qBJ+zjPz2SPMcvMe127RcRO9V9ythMRPwG8jKrf+Z5j2lrT5y9J82QYl7rhXcBrqALi6yPif4+pPMP2wfSGCW1fP+J1/cZWpYf43ph5/aFq2uV2uLg0Is4Ezphyex485XI9zwPeN2Z+/wWBN1F1udg9Ih48qqtKZr6Hgb9SRMTXqCrNO64gYh/go8CTptzmSfs4z89k3MW+814XVNc33P9zfUHmu4BfmrAdPWv9/CVpbuwzLnVAZm5l281/DqIaBWPql0+YP81IE5O6xSxVRBzPtiB+NVV19IepuqnskpmRmUF1M55F6+928YSRS63dW9kWxP8aeC7V6DAPAXbq28ev18us0oghv8C2IP4FqhFYDqeqjvd//r0vVKv03kgqjGFc6o73UgVPgN+acEOT/rs/PmJCu/3zx941siAvqaffAZ6cmX+Qmf+ambcM9Kcf13VlpMzc2At0wx4Di/9j37/H3oRnWhGxG/D/1D9+KDOfm5l/nZlbMvPO+uLNnn3msc6W6X3+XweOrj+vTZl52zw+f0maJ8O41BF1t5Q31D/uz7Z+3sP092k+akLT/fP/dR2b1oQfqqcXZeaNY5Z74hK25c/Z1t3iF+d0189D2HYh6siLHCPiEKq/Bqya3uf/N/UIPaMcuYyNkaRxDONSt/wx226I85tUXRaGuYhqnGYY06UlIvYE/nP949fGjKRSmt71MLuNWiAinsDkLyIzy8ybqPovQzXayTvm0Gz/9T4j95HxX8i6bJrP/6eYMEKNJC2DYVzqkPpP8GfWPz6C6rbzw5b7FvAX9Y/HRsQOoa2+CO4PgH3rp3YYi7xgvTGmn1pXh7dTD3P4wSVuz2+xre/4iyLijyNibBeJiHgw1Sgsw2xm2wWLLxp2B8l6JJHT17m9bdf7/H8yInYYKSUiHkP1uy1JjTOMS93zZ1QXrUE1iscor6Qa6QPgvIh4f0Q8KyKOiIgXAp+iupU6wKXAeQvZ2sX4QD3dDfhkRLw8Ip5SP15DFYz/I9V+LVw9gspz2NY96L8CX4+IP4qIn4+IH4uIH4mIoyPi5yLiXKohJH+wXv57A+19F/ib+sdnAP8QET8bEUdGxIkR8Uf1/G+y/Wg4q6L3+T8SuDQiXhwRR9Xv8+uBz1H937iisS2UpJpDG0odk5lZB46/mrDcdRFxHPC3VH+u/4X6MegzwHOnuJFQMTLz/Ih4H9UoGj/AjlX9e6m+jOwDHL2kbfpmRBwNvBl4MdXIHi9h28WGw9xJNeLHsCEaT6W6Wc2BVDe/OW5g/jVUI6yM/T3oqHdQ3XH0J6j61w/eTOkuqr8aPQP7jUtqmJVxqYMy8yPA5VMs9wWqGwX9KnAx1Wgp91BVU/+OqoL7tMxsyygq98vM/w6cBHya6sY7W4FvUHVPeUpmvq2Bbbo9M0+lGjv8NVTv8ZZ6++6hulvmF6kqu78IPDIzT83Mbw9p69+oguSbga9Q7d8tVH8V+R3g8QM3H1oZmfl94KeBXwH+CbiD6ovNZuAPgSMzc9wdRyVpaWL7EbAkSZIkLYuVcUmSJKkhhnFJkiSpIYZxSZIkqSGGcUmSJKkhhnFJkiSpIYZxSZIkFSsi9o+I90bEdRGxNSK2RMQ5EbHPGtp4ZkS8JSIujIibIyIj4jNTvvY5EfH3EXFjvf5rIuKvI+LJ69+rbbzpjyRJkooUEQcDlwD7AR8BNgFPAk4HToyIYzLzpjFN9JxKdSO0u4CvUd30bdK6dwLeRXVztmuAv6C6c/UjgCdT3evhsjXu0o7rcZxxSZIklSgiPkZ1R93TMvPtfc+/lepOyn+YmS+dop2jgVupwvyjgK8DF2fmU8e85leB/0l1s7gXZ+bdA/MfkJn3rH2vBtZjGJckSVJpIuIgqjvnbgEOzsz7+ubtAXwLCGC/zLxjDe0eyIQwHhF7Av8GfBc4JDO3rm8vJrPPuCRJkkp0XD29oD+IA2TmbcDFwEOouozM23OA3YE/BXaKiBdExG9ExKkR8SPzXJF9xqcUEV8H9qT6diZJkjQPBwK3ZuZjmtqAiPgTYMOCmj+QEdkpM4+c8NrH1tOrRsz/KlUXlsOAC9exbeP8aD29B7gSOKB/ZkR8GPiFzPzerCsyjE9vz53Y5aG777LPQ5veEEmS1A23f/873Mf3m96MDQ96EEccdMh8Y+HVX/s+d93FLGF1r3p6y4j5vef3nmEdo+xXT38N+GfgPwNfBv4j8E7gZ4HbgZNnXZFhfHpbdt9ln4c+5WEvbHo7JElSR1zy7T/n1u/fuKXp7TjokF34i797+FzbfP6zb+TLX/z+pikq4OsV9XQRF0DuXE/vBH46M6+vf748Ip5DVa0/KSJel5n/NsuK7DMuSZKkEvUq33uNmL/nwHLz9J16ellfEAcgM78FfJYqRz9x1hUZxiVJklSir9TTw0bMP7SejupTPo91f3fE/F5Yf/CsKzKMS5IkqUQX1dMT6hvw3K8e2vAYqm4kM994Z4jeBaE/NGJ+7/kts67IMC5JkqTiZOZm4AKqEVlOHZh9FrAb8IH+McYjYkNEzDwyTGb+C9XQiYdHxIv759U/H041Bvo/zbouL+CUJElSqV4GXAKcGxHHUw0zeBRwLFX3lNcNLH9lPY3+JyPiqUAvVO9eTw+NiI29ZTLz5IG2fhH4DPDuiHg+8CWq0VSeDXwPODkz713vjvUYxiVJklSkzNwcEU8E3gCcSBWEvwWcC5yVmTdP2dQhwIsGnttv4LmTB9b9lYg4AjgD+AngGcDNwP8H/E5mXskcGMYlSZJUrMy8BjhlymVjxPMbgY3rXPeLJy44A/uMS5IkSQ0xjEuSJEkNMYxLkiRJDTGMS5IkSQ0xjEuSJEkNMYxLkiRJDTGMS5IkSQ0xjEuSJEkNMYxLkiRJDTGMS5IkSQ0xjEuSJEkNMYxLkiRJDTGMS5IkSQ0xjEuSJEkNMYxLkiRJDTGMS5IkSQ0xjEuSJEkNMYxLkiRJDTGMS5IkSQ0xjEuSJEkNMYxLkiRJDTGMS5IkSQ0xjEuSJEkNMYxLkiRJDTGMS5IkSQ0xjEuSJEkNMYxLkiRJDdml6Q2QpFVx94b97//3rpuubXBLJEmlsDIuSUvQH8SH/SxJWk2GcUlasFHB++4N+xvKJWnFGcYlaYGmCdsGcklaXYZxSVqQtYRsA7kkrSbDuCTN2Xq7nxjIJWn1GMYlaY5mDdT2I5ek1WIYl6Q5mWeINpBL0mowjEvSHCwiPBvIJan7DOOSNINFdysxkEtStxnGJWmdlhWU7UcuSd1lGJekdTAcS5LmwTAuSS3hFwBJqygi9o+I90bEdRGxNSK2RMQ5EbHPGtp4ZkS8JSIujIibIyIj4jNr3I7frl+XEfGMte/JcLvMqyFJWhWGYklajog4GLgE2A/4CLAJeBJwOnBiRByTmTdN0dSpwHOBu4CvAVMH+Xo7jgB+G7gd2H0tr53EyrgktYhfBCStmPOogvhpmfm8zPyNzDwOOBt4LPDGKdt5M/A4qiD902vZgIh4EPBB4P8Cf7mW107DMC5Ja2AYlqTliIiDgBOALcA7B2afAdwBnBQRu01qKzMvzcwvZea969iU3wMeA5wM3LeO149lGJeklvELgaQVcVw9vSAztwvBmXkbcDHwEODJi9qAiDiWqkvMb2bmVYtYh33GJWlKhmBJWrMNEXHFsBmZeeSE1z62no4KwV+lqpwfBly4vs0bLSL2AjYCnwbOnXf7PVbGJamF/GIgaQXsVU9vGTG/9/zeC1r/24F9gVMyMxe0DivjkjQNw6+kLvvufQ/h/FuPmHObnwZu3TRFBXy9op7OPShHxPOBk4BTM/Pqebffz8q4JLWUXxAkdVyv8r3XiPl7Diw3FxHxUOAPgU8AfzDPtocxjEvSBIZeSWrEV+rpYSPmH1pP531h5aOBh1FdQHpf341+EnhRvczH6+deMevK7KYiSS1294b92XXTtU1vhiQtwkX19ISI2Kl/RJWI2AM4BrgTuGzO670J+F8j5j2N6kvA3wPXAV+cdWWGcUkaw6q4JDUjMzdHxAVUI6acSnVBZc9ZwG7AH2bmHb0nI2JD/dpNM6z3GuDFw+ZFxEaqMP7WzPyH9a6jn2FcklrO6rikDnsZcAlwbkQcD1wJHAUcS9U95XUDy19ZT6P/yYh4KtsCdu929ofW4RqAzDx5nhs+LcO4JI1gVVySmlVXx58IvAE4EXg28C2qcb/Pysybp2zqELb19+7Zb+C5k2fb2vUxjEtSB1gdl9RVdbeRU6ZcNkY8v5HqBj6zbsvJzDm0O5qKJA1hVVyStAyGcUnqCL9ASFL7GMYlaYChVpK0LIZxSeoQv0hIUrsYxiVJkqSGGMYlqU8XKstd2AdJWhWGcUmSJKkhhnFJqnWpotylfZGkLjOMS5IkSQ0xjEsS3awkd3GfJKlrDOOSJElSQwzjklZelyvIXd43SeqCTofxiDgpIrJ+vLjp7ZEkSZL6dTaMR8SjgLcDtze9LZLKtQqV41XYR0lqq06G8YgI4H3ATcC7Gt4cSZIkaahOhnHgNOA44BTgjoa3RVKhrBhLkprWuTAeEYcDbwLelpmfanp7JKkEfvGQpDLt0vQGzFNE7AJ8EPgm8Np1tnHFiFkb1rtdkiRJ0jBdq4y/HngCcHJm3tn0xkgq1ypWildxnyWpdJ2pjEfEk6iq4W/JzEvX205mHjmi/SuAI9bbriRJkjSoE5Xxvu4pVwG/3fDmSCqcFWJJUik6EcaB3YHDgMOBu/pu9JPAGfUy766fO6exrZSkhvlFRJLK0pVuKluB/zVi3hFU/cg/A3wFWHcXFkmSJGmeOhHG64s1h97uPiLOpArj78/M9yxzuySVx8pw9R7suunapjdDkkR3uqlIkiRJrWMYl7QyrIpLkkrT+TCemWdmZthFRZK28YuJJJWh82FckiRJKpVhXNJKsBK8I98TSWqeYVySJElqiGFcUudZAZYklcowLkkrzC8qktQsw7gkSZLUEMO4pE6z8juZ75EkNccwLkmSJDXEMC6ps6z4SpJKZxiXJPnFRZIaYhiXJEmSGmIYl9RJVnrXzvdMkpbPMC5JkiQ1xDAuqXOs8EpSd0TE/hHx3oi4LiK2RsSWiDgnIvZZQxvPjIi3RMSFEXFzRGREfGbM8j8YES+PiL+v17c1Im6KiI9HxPPns2eVXebZmCSp3e7esD+7brq26c2QJAAi4mDgEmA/4CPAJuBJwOnAiRFxTGbeNEVTpwLPBe4CvgZMCvIvB34d+DpwEXA9cADwfOAZEXF2Zr5q7Xu0I8O4JEmSSnUeVRA/LTPf3nsyIt4KvBJ4I/DSKdp5M/A6qjD/KKqQPc7lwNMz85P9T0bE4cBlwCsj4k8y84ppd2QUu6lI6hS7qEhSN0TEQcAJwBbgnQOzzwDuAE6KiN0mtZWZl2bmlzLz3mnWnZl/MRjE6+evBD5U//j0adqaxMq4JGk7dlWRNEcbImJo9Tgzj5zw2uPq6QWZed/Aa2+LiIupwvqTgQtn3tLp3VNPvz+PxqyMS+oMq+KS1CmPradXjZj/1Xp62BK2BYCI2BP4WSCBC+bRppVxSdIOrI5Lq+XWex7Ex6/fMOc2/wm4ddMUFfBR9qqnt4yY33t+73W2vyYREcB7gEcA59VdVmZmZVySJEltFPU0l7S+twAvBD4NzGUkFTCMS+oIu6jMn++ppIb1Kt97jZi/58ByCxMRv081esungGdn5tZ5tW03FUmSJJXoK/V0VJ/wQ+vpqD7lcxERZwOvoBpv/Kcy83vzbN/KuKTWs4IrSZ10UT09ISK2y6wRsQdwDHAn1bjfcxeVd1IF8Y8DPznvIA6GcUnSGH7RkdSUzNxMNWLJgVR30Ox3FrAb8IHMvKP3ZERsiIiZr0StL9b8I+BlwN8Dz8nMO2dtdxi7qUiSJKlULwMuAc6NiOOBK4GjgGOpuqe8bmD53ggn0f9kRDwVeHH94+719NCI2NhbJjNP7nvJ6+vl7wQ+D/xGlc+38/nM/Ks179EAw7ikVrNyu3gOcyipKZm5OSKeCLwBOBF4NvAt4FzgrMy8ecqmDgFeNPDcfgPPndz378fU0wcDvzmizfcDhnFJkiR1V2ZeA5wy5bI7lK/r5zcCG9ewzpPZPpwvjH3GJbWWVXFJUtsZxiVpiW45+IHccvADm96MNfOLjyQthmFcUiu1LRwOhvA2BnJJ0vwZxiVpwUYF77ZVydv2BUiS2sALOCVpQaYN2r3l9to8t7srS5JawjAuqXVKr9Cut9ptKJek1WM3FUmak3l1Oym560rpX4QkqW2sjEtqlRLD4CLCs1VySVoNVsYlaQaLrmKXeJFniV+IJKmtDOOStA7LDsmlBXJJ0nwYxiW1xqpXZEsK5Kv+WUjSvBjGJWmNmgzFJQVySdLsDOOSWqGUSqxhWJI0T4ZxSWqZUr4QlPIFSZLazDAuSVMqJQRDWdsiSVo/w7ik4pVQgTX8DlfCZyNJbWYYl6SW8guCJLWfYVxS0UqovJYcekveNknSZIZxSRqjDWG36W0s4QuTJLWVYVySRmg65EqSus8wLqlYVlyn1/QXBz8rSVofw7gkDdF0uF2PNm6zJK06w7ikIjVZaTXUSpKWxTAuSR3S5BcJu6pI0toZxiWpTxeq4l3YB0laFYZxScVpqsJqiJ2d1XFJWhvDuCR1kF8sJKkdDOOSimJVfH66uE+S1DWGcUkrr8uhtYl9s6uKJE3PMC6pGE2EuC4HcUlS+QzjktRxVsclqVyGcUlFsCouSVpFhnFJWgFWxyWpTIZxSY2zKi5JWlWGcUlaEVbHJak8hnFJK8equCSpFIZxSY2ycrpcVsclqSyGcUkrxaq4JKkkhnFJjbFi2gyr45JUDsO4pJVhVXwb3wtJbRER+0fEeyPiuojYGhFbIuKciNhnDW08MyLeEhEXRsTNEZER8ZkpXvcfI+LPIuKGiLgrIr4SEWdFxINn26ttdplXQ5K0FlZKJUmTRMTBwCXAfsBHgE3Ak4DTgRMj4pjMvGmKpk4FngvcBXwNmBjkI+Io4BPAA4DzgWuA44DXA8dHxPGZuXXNOzXAyriklWAleEfLfk/8AiZpHc6jCuKnZebzMvM3MvM44GzgscAbp2znzcDjgN2Bn560cETsDLwPeAjwgsz8+cz8deAo4MPAMcAr17ozwxjGJS2doUySNElEHAScAGwB3jkw+wzgDuCkiNhtUluZeWlmfikz751y9T8OHA58KjP/uq+d+4Bfq398aUTElO2NZBiX1HlWxUezOi6pYMfV0wvqEHy/zLwNuJiqcv3kBa77o4MzMvNq4CrgAOCgWVdkn3FJnWYQl6RGbYiIK4bNyMwjJ7z2sfX0qhHzv0pVOT8MuHB9mzfTug+rH5tnWZFhXNJSWRktzy0HP5C9Ns98DdLU7t6wP7tuunZp65M02d13P4At1z587m3OaK96esuI+b3n9551RU2u2zAuqbOsiktS4zZNUQFfr15/7VxQ+0tZt33GJS2NVfFy2XdcUoF61ee9Rszfc2C5Vq67U2E8It5cD+Z+TUTcWQ/q/s8RcUZE7Nv09klaHqvia+d7JqkwX6mnh42Yf2g9HdWvuxXr7lQYpxrvcTfg48DbgD8Bvg+cCXwhIh7V3KZJq81KqCRpjS6qpydExHaZNSL2oBrr+07gsgWs+xP19MTBGfWQi4cB3wCunnVFXQvje2bmkzPzv9eDwr88M38U+F3gB4DfbHj7JC2BFd71W+Z75xc0SeNk5mbgAuBAqjto9juLqgD7gcy8o/dkRGyIiA1zWP0ngSuBp0XEc/ra34nqBkIA78rMmfuMd+oCzsy8a8SsPwNey7Y/KUhaIkOXJGmdXgZcApwbEcdTBeSjgGOpuoi8bmD5K+vpdjfjiYinAi+uf9y9nh4aERt7y2TmyX3/vjciTqGqkJ8fEecD3wSOB55INcb52TPuG9CxMD5G77anX2h0KyQtnFXx2S1zqEOHOZQ0TmZujognAm+g6jLybOBbwLnAWZl585RNHQK8aOC5/QaeO3lg3Z+NiB+lqsKfAOxB1TXlDcCbMnMuB8pOhvGIeA3Vt569qL69PJUqiL9pitcOHZgemMefPCQtkEFckronM68BTply2aG3p8/MjcDGdaz7y8AL1/q6tehkGAdeAzyi7+ePAidn5o0NbY+0suyisr3bDqjOE3t8o4lhcadndVySlqOTYTwzHwkQEY8AnkJVEf/niPipzPzchNcOHZi+rpgfMe9tlTQfJVfFewF82HOlh3JJ0mJ1bTSV7WTmv2fmX1L189kX+EDDmyStlFWvit92QAwN4sOWmbRcExxZRZIWr5OV8UGZ+Y2I+DLw+Ih4WNEgu88AACAASURBVGZ+u+ltkjQ/pVXF1xusrZZL0urpdGV8wA/U03sb3QppRaxapXOeFe6SKuVWxyVpsTpTGa8HeP9uZl4/8PxOwO9QDV9zSWZ+p4ntk1bJMkNV01XxRYbm/ratlktSN3UmjFONPfn7EfEpYDNwE9WIKj8OHARcD7ykuc2T1CXLrlzfdkA0FsgdWUWSFqdLYfwfgD8CjgF+BNgbuIPq7kwfBM5dw8DwktZpFariTXUhWZVALkmrpDNhPDO/CJza9HZI0iI1GciXxeq4pFWyShdwSlowq+Ld3oam++dLUhcZxiVpSiUE8Z6StmURHFlF0qowjEuai65XxUsMv01sk9VxSZovw7gkTVBiEO8pedtmZXVc0iowjEuaWZer4m0Iu8veRqvjkjQ/hnFJrWEQXz1WxyV1nWFc0ky6GpbaFsStjktSOxnGJbXCMsNf24J4T1u3e5KufuGTJDCMS5qBIak8ywzkVsclaXaGcUnFsyq+Nl3Yh0F+8ZPUVYZxSevSxXDUpRC7rH2xOi5JszGMSyrassJel4J4T9f2qYtfACXJMC5pzboWiroWWvstY9+WWR3v2u+eJBnGJRVrGSGvy0G8p2uBXJK6xDAuaU26VJlchSDeRV36HZQkw7ikIllpnS+r45JUJsO4pKl1qSK5ilXxLu1zl34XJa02w7ik4lhhbS8/O0laG8O4pKksqxLpRZuL1aV9tzouqQsM45JWSpfC6Hot+j2wOi5J0zOMS5qoK1Vxg3j3WB2X1HaGcUljGXa6yeq4JJXBMC6pCFbFl68rgdwvjJLazDAuaaSuhByDePd15XdV0uoxjEsaapnhxi4NzelKdVyS2sowLqlRdk9pXlfeI6vjktrIMC5pB10JNV0JmW1ndVySRjOMS2qMIa0cXfni0pUvkpJWh2Fc0na6Ema6Ei6XaZHv2TK/eHXld1jSajCMS7pfVy7aNIiXyb+ESNKO5hbGI+KPI+Jl82pPUncZysrVlS8yVsel7oiI/SPivRFxXURsjYgtEXFOROyzxnYeWr9uS93OdXW7Iw8YEfGTEXFBRFwbEXdGxNUR8ecRcfTse1aZZ2X854FzI2LXObYpaUm6El66Eiab1JXuKpLaLyIOBq4ATgEuB84GrgZOBy6NiH2nbGdf4NL6dZvrdi6v270iIg4a8po3A/8HOAL4KPA24HPAc4GLI+K/zbRztV3m0Ugfz4KSxlrl7ilbH333/f9+4DfLrlvcdkCwxzey6c2Yyd0b9mfXTdc2vRmSZnMesB9wWma+vfdkRLwVeCXwRuClU7Tzu8BhwNmZ+aq+dk6jCtnnASf2Pf9I4DXAvwM/nJk39M07FvgE8Abgj9e9ZzX7jEvqRFW81CC+9dF33/+Y5vlV4MWckqZRV6tPALYA7xyYfQZwB3BSROw2oZ3dgJPq5c8YmP2Ouv1nDVTHD6DKyZ/tD+IAmXkRcBvw8DXszkiNhfGoPCMi5rIjktanKxdtlmYtQbvUYG53FUkNO66eXpCZ9/XPyMzbgIuBhwBPntDO0cCDgYvr1/W3cx9wQf3jsX2zvgrcDTwpIh7W/5qIeBqwB/AP0+/KaPPupgIw7d81d6Xa+fsWtB2SVkQpVfF5hOk2dWVpC7urSI3aEBFXDJuRmUdOeO1j6+lVI+Z/lapyfhhw4YztULfT27abI+LXgbcCX46IvwJuAg4GngN8HPjlCds/lUWE4M9HxBeALwD/CvxrZn5jyHKPrKf3DZknaQmsis9ukdXsEoL5IvuO33LwA9lr89aFtC2pE/aqp7eMmN97fu9FtJOZ50TEFuC9wEv6Zn0N2DjYfWW9FhHGDwc2AP+590RE3AZ8sX58larPzgvq2dctYBskFaSLF20uu0tJb31NhHIv5pS6L+6OuR9f4u4A2DRFBXzdq6insx6ghrYTEb9GdeHnuVR9y6+nyri/B/xJRDw+M39txnUvJIyfTRXIfxj4gfq5PYGnUPXZGfShBWyDpAm6cGHbqgTxwXV3qfvKMqvjBnKpdXoV671GzN9zYLm5tRMRTwfeDPxl/+grwOci4meoury8OiLelZlXT1j/WIsI47+ZmXcD1IOx/3Df44eAg4DdgRuB84HfWsA2SBrD7inrV8JFlk0E8i5UxyW1zlfq6WEj5h9aT0f1BZ+lnZ+qpxcNLpyZ34uIy4GfAZ5ANe75ui30wsnM/A7wyfohSXPTRFW8hCDe06VAbnVc0gi9IHxCROzUP6JKROwBHAPcCVw2oZ3L6uWOiYg9+kdUiYidqC4C7V8fQK+SNGrUv97zM58YHGdcWjFdqIqvehDvKXGb1qtrf0GRNLvM3Ew18t6BwKkDs88CdgM+kJl39J6MiA0RsWGgnduBD9bLnznQzq/U7X9soLvJp+vpL0XED/a/ICJ+guqLwF3AJWvdr0EOKShpIboSrkoPvMuukHehu4rVcalVXkYVeM+NiOOBK4GjqMYEvwp43cDyV9bTwarNa4GnA6+KiMcDl1Nd4/hc4AZ2DPvnU40j/gzgyoj4S6oLOA+n6sISwG9k5k0z7t9cK+O/DvwVDlUoFcuLNtem9CDes+ztXNRn4J05JQ2qq+NPBDZShfBXU431fS5w9LRhuF7u6Pp1h9TtHAW8DziyXk//8vcBzwZeCXyZqn/4q6luMPR3wLMy820z7h4wx8p4Zv7+vNqSNH92T1mbtgTxnq6MsmL/cUmDMvMa4JQplx15osjMm4HT68c0bd0DnFM/FsY+49IK6EIQX6a2BfGerY++e2nbXspdTyWp7QzjklphWeGvrUG8X9sDud1VJK0Sw7jUcV2oihvE167t+2Igl7QqDONShxkyptf28DrMMvbJ7iqSNBvDuKS5aHNVvItBvKfNgdzquKRVYBiXOsruKdPpchDvWYV9nAcDuaQmGMalDjJUTGeVQuqi97UL1XFJaoJhXNJM2loVX6Ug3mMgn8wvspKWzTAudYzdUyZbxSCu6RnIJS2TYVzqEEOEJrE6Ph3/L0laFsO41BHLDg9WxdvLQC5J5TCMS1ozg3j7+V5MZnVc0jIYxqUOMDRMZvjc0SLfk65Ux/2/JWnRDONSy9k9ZTKDeDO6EsglaZEM45KmZhDvHt+fyayOS1okw7jUYoYEzYPdVSbz/5qkRTGMSy1l95TJrPpOz0AuSc0wjEstZBCfzCC+dr5n41kdl7QIhnFJnWOoXL9FvXddqY4byCXNm2Fcahmr4uMZxGdnIB/PQC5pngzjUosYxMcziJfPQC5J2zOMSy1hEB/PID5fvp+TGcglzYNhXGoBT/pqgt1VJGnxdml6AySN10QQtyo+Xwfuf+PY+VuuffiStmTttj76bh74zV3n3u5tBwR7fCPn3u4tBz+QvTZvnXu7o9y9YX923XTt0tYnqXsM45K2YxCf3aTwPWn5ksN5GxjIJbWJYVwqWFf6iS9KSUF8rQF82rZKCOZtq443wUAuab0M41KhutRPfJE392nKPMP3tOtpMpi3LZAvuzoOBnJJ6+MFnFKB7Cc+WZNV8WUF8WHrbWrd4AWd0+jSl2hJy2EYlwpjEJ9sFYN4KdtQUtegabSt65Wk1WMYl1Zc28JKU2Gw6ar0oJK2ZR661JXJ6riktTCMSwXp0km8S+Gq1ODb1BcEu6tM1qX/y5IWqzNhPCL2jYgXR8RfRsTXIuLOiLglIj4TEb8YEZ3ZV3WT3VMma6IqXmoQ72cgn8xALqlUXQqoLwTeDRwFfBY4B/gw8DjgPcCfRUR3SnXqFIP4ZMsO4qV1S5mkS4F8UQzkkkrUpaENrwKeA/xtZt7XezIiXgtcDvws8HyqgC4Vo0tBfFGaCOJt1NvuEsYmn0WXxh+XpEk6UxnPzE9k5t/0B/H6+euBd9U/Pn3pGyaN0bWq2SKq4gbxtVvmPthdZbKu/T+XNF9dqoyPc089/f6kBSPiihGzNsxvc6TmTtBt656yLF0I4f0O3P/GpVXIvSHQZN4QSNIonamMjxIRuwC/UP/40Sa3ReoxiE9nWVXxrgXxnmX2e7f/+GRWyCUN0/kwDryJ6iLOv8vMj01aODOPHPYANi18S6UFsp/4cF0N4v3aHMgX+RcXA7mkEnQ6jEfEacCrqYL0SQ1vjgR074LNNvcTX4Ug3mMgL4eBXFK/zobxiDgVeBvwZeDYzLy54U2SDOJTMIgvTpv3uUsXdIKBXNI2nQzjEfEK4B3AF6mC+PUNb5LUuZNvm6uVbQ6ls1rGvtt/fDpdOyZIixIR+0fEeyPiuojYGhFbIuKciNhnje08tH7dlrqd6+p2x/5njIgfi4gPR8S36td9KyIuiIhnz7Znlc6F8Yj4deBs4PNUQfyGhjdJ6twFm4uyjBC3ykG8p62BvGv9x8FALk0SEQcDVwCnUN035mzgauB04NKI2HfKdvYFLq1ft7lu5/K63Ssi4qARr/st4FPA06gGAnkL8DfAPsxpyOxODW0YEb8NvIHqQzvBrikqQReDeFu7pxjEt1nG0IeLGPJwkTcEamLIQ3DYQ2mC84D9gNMy8+29JyPircArgTcCL52ind8FDgPOzsxX9bVzGlW35vOAE/tfEBEvBH4H+Afg+Zl528D8B6xnhwZ1pjIeES+iCuL3Ap8GTouIMwceJze6kVo5BvHpGMSb0db3xAq5tBrqavUJwBbgnQOzzwDuAE6KiN0mtLMb1UAed9Sv6/eOuv1n9VfHI2In4M3A94CfHwziAJl5z+Bz69Glyvhj6unOwCtGLPNJYONStkYrzyA+HYN4sxZdIW/bDYGaZIVc2sFx9fSCIXdYvy0iLqYK608GLhzTztHAg+t2tgvVmXlfRFwA/BJwLFUXGICnUGXL84HvRMRPUg2VfRdweWZeOtOe9elMGM/MM4EzG94MCTCIl8QgPllbA/miNNVdBQzk6qQNo+5uXt/HZZzH1tOrRsz/KlUYP4zxYXyadqjb6fnRevrvwOeA/9T/goj4FPCCzJz5JNOZMC6VootBfFEWXRU3iE+vjYG8i/3HwUCuZux8N3P//7Tz7If4verpLSPm957fewHt7FdPXwp8HXgG8FngAKqLOJ8F/DlzuIizM33GpRJ0td9nG7unGMTXbtHvmSOsTK+rxxKtpE1j7m4+q94BYNZvEcPa2blv3gsy88LMvD0zvwT8DHAt8OMRcfSM6zaMS/PS5Mmzbd1TDOLlauN7ZyCXOqtXsd5rxPw9B5abZzvfqadXZ+a/9C+cmXcCH6t/fNKEdU9kGJfmwCA+PYN4+Rb5Hi7q8zeQS530lXp62Ij5h9bTUX3BZ2mn95rvjnhNL6w/eMK6JzKMSzMyiJfDID4/BvLtGcilRlxUT0+ohxq8X0TsARwD3AlcNqGdy+rljqlf19/OTlQXgfavD6ob/XwfODQihl3w8rh6umXCuicyjEszMIivzSKr4gbx+WtjIF8kA7m0XJm5GbgAOBA4dWD2WcBuwAcy847ekxGxISI2DLRzO/DBevkzB9r5lbr9j2Xm1X2v+TbwIaquLa/vf0FEPJPqAs5bqO7KORNHU5HWySC+NgZxDWrbCCvgKCtSA14GXAKcGxHHA1cCR1GNCX4V8LqB5a+sp4Mns9dSjXzyqoh4PHA5cDjwXOAGdgz7AK+q1/W6iHha/ZoDqC7gvBd4SWaO6sYyNSvj0jp0NYgvikG8vRxhZUdWyKXlqavjT6S6aeNRwKuBg4FzgaMz86Yp27mJ6uY/5wKH1O0cBbwPOLJez+BrbqiXORt4FHAa1Y2I/hb4scz881n2rcfKuLQGTZ8IFx0C2njBphZv0WOQL0IX79DZY4VcqyYzrwFOmXLZkSeyzLwZOL1+TLvum6kq5K+a9jVrZWVcmpJBfO0cOaU72th/vKsXdELzxyNJ82MYl6bQ9ImvjUF80doaxJ/5yE0885Gbmt6MdTGQb89ALmke7KYiTdD0Ca+tQdx+4owN3YPzPn79hhFLlmWRXVYWcUEnLLbLSpMXdIJdVqQuMIxLYxjE12dVg/gsFe9hry01oBvIt1dCIAcM5VJLGcalEQzi67NKQXzR3U1Krp63MZAvUtOBHKySS21ln3FpCIP4+qxKEG+q33dp/c3b1oe8y0Me9jR97JK0doZxaUDTJ7MSTujrsUpBvGklbENPSZ/NNAzkkkpjGJf6NH0SW8aJvG1DGJYS9kqrSpe0PYv6jNo4wgoYyCWtjWFcojpxNX3yMojvqIQgXlLoHabkbZsHA/n6NX1MkzQdw7hWXgknrLYG8UUqJYi3QQlfGNrWfxxWJ5CXcIyTNJphXCuthJNUm4P4okJS00G8hHC7Hk1vt4F8RyUEcijjWCdpOMO4VlYJJyeD+I5KCOJtZyBfm2UE8hJCeQnHPEk7Moxr5ZTyZ1uDeFmarirPW5P7YyAfzkAuaRjDuFZKKScig/hwTVXFuxTCBxnIp7dKgbyUY6Ekw7hWSCknH4P4cE0E8a5Vw0cxkE9vVQI5lHNMlFadYVwroZSTjkF8uKaC+CoxkE/PQC5pmQzj6rSS/hxrEB/OIL48BvLpGcglLYthXJ1VyglmWSMpGMSns6pBvKeprjlNj5KzHqsWyEs5ZkqrxjCuTirlpLKsE61BfDqrHsT7dSmQL/L3dJUCOZRz7JRWiWFcnVJSdccgPppBvAwG8uksK5CXEspLOYZKq8Iwrs4o6QRiEB/NIF6WLr03bQ/kUE6VvKTChtR1hnF1QkknDYP4aMsO4qsydOGslv0etfGCTli9QA5lHVulrjKMq9VKq94YxEdrIohregby6axqIC/pOCt1jWFcrVXSyWGZI6YYxCcziK/Psv+SYCAfr6RADmUdc6UuMYyrdUqr0rS9Gg4GcW3PQD7ZMgN5SaG8pGOv1BWGcbVKaScCg/h4BvH2MpBPtqxADmVVyUsriEhtZxhXa5R28DeIl8UgPn8G8skW2XVsUEmBHMo7JkttZRhX8Uqrwizzz8ZtDuLLrIobxBfHQD6dVQ7kJR2fpTYyjKtYJR7kl3kiNIhPxyC+eAby6axqP3KwSi7NwjCuIpV4YDeIT8cg3k0G8umsaj9yKLOAIrWBYVxFKfFgvuwqlEF8Ogbx5TOQT2eVAzmUWUyRSmYYVzFKPIAvO4QbxKdjEG9OlwJ5l0ZaKS2Ul1hYkUplGFfjSj1od6UaDgZxzVdXAjl0Z+hDKLdKXuLxXSqJYVyNKfUg3aVuKWAQ12IYyKezzKEPocxADmX+5VMqhWFcjSj1wLzsE5lBfHoG8fIYyKe36t1WoNwCjNQ0w7iWquSDsUF8bQziAgP5WthtpVLqOUBqimFcS1F6CO/ShZpgENdyGcin10QgLzGUl3xOkJbNMK6FK/mA27VqOBjE1QwD+fSW3Y8cyq6Sl3yOkJbBMK6FKfkg20S1aBnVcIO4mmQgXxur5NuUfL6QFs0wrrkr/aDaRAhve7cUMIhrOgbytVl2IIdyq+RQ/vlDzYiI/SPivRFxXURsjYgtEXFOROyzxnYeWr9uS93OdXW7U/3SRcRJEZH148Xr25sdGcY1N6UfRLtYDQeDuMrTtUBut5XlK/lcouWKiIOBK4BTgMuBs4GrgdOBSyNi3ynb2Re4tH7d5rqdy+t2r4iIgya8/lHA24Hb17cno+0y7wa1etpw0GzixGMQX7suB/EX7Pm5oc+ff+sRS96S5XjmIzfx8es3LGVdB+5/I1uuffhC17H10XfzwG/uutB13HZAsMc3cqHr6Nc7Lu61eevS1rkWvXPLrpuubXhL1LDzgP2A0zLz7b0nI+KtwCuBNwIvnaKd3wUOA87OzFf1tXMa8LZ6PScOe2FEBPA+4CbgL4DXrGtPRjCMayalB/Gmqj8G8bXrWhAfFb4nLdelcG4gX7tlB3KojpOlBnIwlK+yulp9ArAFeOfA7DOAXwJOiohXZ+YdY9rZDTgJuKN+Xb93UIX6Z0XEQZl59ZAmTgOOA55eT+fKMK51KT2EQ3er4dCtEVOg/UF82uC93rbaHNAN5GvXVCCHcqvkYChfUb3ge0Fm3tc/IzNvi4iLqcL6k4ELx7RzNPDgup3bBtq5LyIuoAr2x1J1gblfRBwOvAl4W2Z+KiIM42qWIXw0q+Hr09YgPs8AvpZ1tTGYG8jXrnc8MZTvyFDeOhsi4ophMzLzyAmvfWw9vWrE/K9ShfHDGB/Gp2mHup37RcQuwAeBbwKvnbCt62YY11QM4aN1pRoOBvFpLDOEj1q/gXy8ZQVyoJNVcii/6woYyudt5633zf0z33nrfZMXGm+venrLiPm95/deUDuvB54APDUz75ywjnUzjGusNoRwMIjPg0F8sqaDeE9vO9oWypcdyAGr5DNoQ5UcDOUtsGmKCvh69U7Cs/7n2KGdiHgSVTX8LZl56Yztj+XQhhqq9GEKe5q6icUyhyIziDfvBXt+rpgg3q/EbZpk2Z/9Mn63l/F/FJoZkxzKvllQv7act7QmvYr1XiPm7zmw3Fza6euechXw25M3czaGcW2nLQezJk8OywzhBvHmlR54S/2iMI6BfP2aGJO8pw2BHNpzHtNUvlJPDxsx/9B6Oqov+Hrb2b1e9nDgrr4b/STbRmN5d/3cORPWPZHdVAS0pzsKNHtC6FI1HAzi47Qt4LatL/kyu6xAt/qRQ7N9yaH8ritg95WOuKienhARO/WPqBIRewDHAHcCl01o57J6uWMiYo/+EVUiYieqi0D717cV+F8j2jqCqh/5Z6hC/sxdWAzjK84QPp1lVqIM4s1rWxDvaVtf8i4GclhOP3Jori85GMq1HJm5uR528ATgVKo7YPacBewG/GH/GOMRsaF+7aa+dm6PiA9SDV94JvDqvnZ+BTgQ+FhvjPH6Ys2ht7uPiDOpwvj7M/M9s+1hxTC+ogzh0zOIz64tQbytIXxQm6rkBvLZNVUlh3aMutJjKG+tlwGXAOdGxPHAlcBRVGOCXwW8bmD5K+vp4Mn7tVQ37XlVRDweuJyqG8pzgRuown4j7DO+Qnr96Azi01n2RZoG8WZ1JYj3tKkveRf7kMPyvlxD833Jmy6arEXbzoOrLjM3A08ENlKF8FcDBwPnAkdn5k1TtnMT1c1/zgUOqds5iuo290fW62mElfEV0MaDTtMHdqvh82EQb15bquRdrpDDcvqRQ/NVcmhH1xWwUt4mmXkNcMqUy448gWfmzcDp9WO923ImVVeXubEy3mFt/PbfdIVl2dUlg3jzuhzEe9qyj01UyK2Sz1/Tx/G1auNfjdUthvEOauNBpemDdxMh3CDevLaE1Hloy7428buzzEC+7FDepKaP6+vRxvOn2s8w3hFt/WZfwsF62SesZZ6MDeKjtSWczlNb9rnLgRxWq0oOZRzn16qN51O1l2G85dp6wCjh4NzlajgYxMdpSyhdhLbsu4F8vgzl69PWQpfaxTDeUm09OJRwMG7ipLTsEG4QH60tYXSR2vIerEIgbyKUN62E88B6tPW8q/IZxlukzd/QSzn4NnEi6no1HAziWpyuB3JYzSo5lHNeWKu2nodVLoc2bIE2/6cv5UDb9RAOBvFJDOLba8uQh7D8YQ9h2/+nZQx/CMsfAhGavYNnv7YNidjTf252eETNwsp4wdr87buUikdTFSCDeFkM4sO16X1p6net61VyKKPrCpRTvFmPNp+v1TzDeGHa3BUFygnh0Fw13CCuNmlTIG/KqgTyEkJ5SeeQ9Wjz+VvNsZtKAbrwH7ekg2dTJ5RVCeHQriBu2JysLV1Wmuiu0rOsO3b2NNFtBcrrugLt674CdmHR2hjGG2QIn69VCeFgEJ+WQXx6BvLJlt2PHAzl0N4+5T0Gc01iGF+yLgRwMIT3GMTLZRDvriYDOSy/Sg7VsWbZgRwM5fNmMNcwneozHhEviIi3R8SnI+LWiMiI+OOmt6vt/cB7en35SgniTfZxbKJvOBjEtVht+gLT9O9jE/8XmzruQDkXeUL7+5X39HJBPugBTW+KGta1yvhvAT8C3A5cCzRWOml78O5X2kFv1Srh0GwIh+aDz1q1KVSWpi3dVWA1K+Rglbyn7f3KpZ6uhfFXUoXwrwE/Dly0zJV3KYBDeSEcDOJNMIirZCUEclhuP3Jori85lBfKoRtdWLS6OtVNJTMvysyvZuZSjxBd6ILSr8Q/Aa5ilxQwiKsZbftCU8LvaVP/V5s6NkE5wyH2K/H8JU3Stcr40nQpfEOZVXBovp/iqoZwKCPgrFXbQmTJ2tRdBZqvkEOz3VagmSo5WCmXZmUYHxARV4yYtSEf9ABD+JKsaggHg/h6GcRVSiCH5XdbAUP5MPYrVxt0qpuKplPaqCj9mv6zZ5NdUsAgrrK08QtOKb+/Tf5fbvIYBs0fx0cp9bwnWRkfkJlHDnu+rpi352+2Q5R8ECrhwN30CayEIC51QQkVcmiu2wo0XyWHMivlYLVc5TGMd5wBfDJD+DalVBWlWRnIK4by8QzmKoFhvKMM4ZM1HcLBIC4tUkmBHJrpR95jKJ/Miz7VFMN4h5QcwMEQ3q+kEA7tD+Jt7Nus5SglkEPzVXJo7oZB/frPBSUGc6vlWjYv4Gy5ki/G7CnpYh6D+I7aHsS1WF34olPS7/iB+9/Y+DGg6QvV+5V0fhimDedYtV+nKuMR8TzgefWPj6ynR0fExvrf387M1yx9w+asLQeFkg6wpZx4mj4JDyoppEiLVFKFHMqpkkOzXVd6Su/CAnZj0eJ0KowDjwdeNPDcQfUD4BtAa8N4G0J4SQEcDOHjGMS1agzkwxnK18ZuLJq3ToXxzDwTOLPhzZirNgRwMISPYxCXylFiIIdmL+7sKTGUg8Fc3depMN4VbQngYAgfp8QQDgZxqUSlVMmhjIs8+7WhWg4Gc62fYbwQBvDZlBTCwSAulay06niPVfLx2hLKwWCutXE0lQa17SrtEq96L2lUgJ5Sg3iXnX9rq2+OW7SuvrclfzEt6RhS4jG2dy4q7Xw0StvO9Vo+K+NL1rb/jCUf7Eo7QZR0Ah2m1TABIAAAGHxJREFU5PAhNaHUCjmU1W0FyqyUQ7uq5WDFXMNZGV+CNn4rLrnqUGKlxiDevK5WcJu0Cu9pyf83ShiTfFCJx19oX7UctmWDex9oFFt1VsYXoE2hu1/pB7ESTwClnSiHKTlsqFyrEMTborQqOZRbKYf2jMQi9fh1bE7aWP3uKb2aUGolpg1BfNUYILUebfjCWmKVHMo9PveUfn6TwMr4TNoYvHtKPziVfHAv8YQ4ShtCxrydf+sRnbiFe5NW8UtNyf3H+5VYJYeyK+VgtVxlM4yvkQF8sUoO4WAQbwsD+fqtYhDvaVMghzKGQBxUeigHg7nKYzeVNWjjRRZtuail9D91QruCuFY7VK6X71m7lHxM6h3TSz+ut+H8qO6zMt5RbTm4lH6ghrJPeKOsclW8X3+4tFI+nAG83UrtttLPark0nmG8Q9oSwKEdIRzaGcQ1XC90GsoN4OO0patKv5K7rfRrQygHg7mWzzDecm0K4NCeEA4G8a5a1Wq5Abz72lAlh/aEcjCYlyIi9gfeAJwI7At8C/gr4KzM/M4a2nko8HrgecB/AG4CPgq8PjOvHVh2X+BngJ8E/hPwg8DdwL8C7wPel5n3zbZnFcN4C7UtgIMhfJnsojK9wYDapXBu+F5NbQnk0K5QDgbzpkTEwcAlwH7AR4BNwJOA04ETI+KYzLxpinb2rds5DPgE8KfABuAU4Ccj4ujMvLrvJS8E/oAq+F8EfBN4BPB84D3AT0TECzNz5l8Gw3hLGMCXo+1BXLMZFWBLD+kG7/lqY1eVfm3pttLTtlAOBvMlO48qiJ+WmW/vPRkRbwVeCbwReOkU7fwuVRA/OzNf1dfOacDb6vWc2Lf8VcBzgL/tr4BHxGuBy4GfpQrmH17fbm1jGC9YGwM4GMLVPaWEdEO31qKtoRwM5qpExEHACcAW4J0Ds88Afgk4KSJenZl3jGlnN+Ak4I76df3eQRXqnxURB/Wq45n5iWFtZeb1EfEuqi8BT8cw3i1tDd/QzgDe06UgbheV5TIcqw3a1HWlp43VcjCYL8Bx9fSCwf7ZmXlbRFxMFdafDFw4pp2jgQfX7dw20M59EXEBVbA/Frh6yOsH3VNPvz/FshMZxhvW5gAOhnBJ3dT2riqD2lYl72lrtRx2PL+vcDjfEBFXDJuRmUdOeO1j6+lVI+Z/lSqMH8b4MD5NO9TtjBURuwC/UP/40UnLT8Mw3oC2B3BodwgHg7ik1dTGKnlPW6vlPVbN12WvenrLiPm95/deUjsAbwIeB/xdZn5siuUnMowvSRcCOLQ/hEO3g/jHr99gVxVpDrpUFR/U5kAO7Q/lUGYmiLvuYddN105ecI1tApumqICvexX1dNZvN1O1U1/s+WqqEV1OmnGd9zOML0iJ/9HWqwsBvKfLQbzHQC7NpstBvKftgRy6EcoB7m335i9ar2K914j5ew4st7B2IuJUqlFXvgwcn5k3T1jn1Azjc9SlAA7dCuGwGkFc0mxWIYj3dCGQQ7v7lWuir9TTUX25D62no/qCz6WdiHgFcDbwRaogfsOE9a3JTvNsbNXcdkBs9+iKrY++2yDecqsUKKR5WcX/N107Nnbx/LXiLqqnJ0TEdpk1IvYAjgHuBC6b0M5l9XLH1K/rb2cnqotA+9fXP//XqYL454Fj5x3EwTC+Zl0M37DtANbFg1jXTjbTWsVgIWntuniM7Or5bNVk5mbgAuBA4NSB2WcBuwEf6B9jPCI2RMR2J8DMvB34YL38mQPt/Erd/scG7sBJRPw21QWbV1BVxL892x4NZzeVNehiv64uH6y6eIJZK/uPS9NZ9S+vXemyMsguLJ3wMqrb2J8bEccDVwJHUY0JfhXwuoHlr6yng1XT11LdpOdVEfF4qrtoHg48F7iBgbAfES8C3gDcC3waOC1ih0LslszcuM79up9hfAV1OYD3GMS3MZBL4616EO9p61jk0zKYt1Nmbo6IJ1IF4xOBZwPfAs4Fzpr2QsrMvCkijqa6A+fzgB8DbgLeB7w+MweHknlMPd0ZeMWIZj8JbJx+b4YzjK+IVQjgPQbxHRnIpeEM4jvqapW8n8G8XTLzGuCUKZcd2Y+4Du6n149J7ZzJjl1aFsI+4x23av3mDOKjGTqk7fl/QrB650mVx8p4R3lg0TC98GGVXKvKAK5RujJuudrHMN4xhnBNY1ggMaCriwzf67MKXVVGMZRr2QzjHWEI16wGQ4vhXG1k+Na8GMq1LIbxFjOAa5EM52oDw7cWzYs9tWiG8ZYxgKspdm1RkwzdKoHVci2CYbwlDOEq0aiAZEjXehm61QZWyzVPhvGCGcDVVpMClWF9tRm41SVWyzUrw3iBDOHqumnCmIG9nQzaWlVWy7VehvGCGMKlbdYT6gzw82WwltbHarnWwjBeAEP4/PTGxfVOnKvJ8CjNblXHF///27v7YEmq8o7j34dVEDeAgFoYQVc2LKsmgq8IVJSXEjFUkCQarYooJJYm8YUkWhXLJAImmkqVJWrUmIQo4hsm5IWKhYLlsmIkxkAUwQAawqIG5FVgF2RX2Cd/9Bl3dvbeO3Nnuqd75n4/VVO9091z+pyzfWd+03O6uwlbn7TNQK6hDOMtM4g3w1AuSaMzgDfHo+QaxjDeIoN48wY/YAznkmT4boNHybUYw3hLDOLtMJxLWqkM4O0zkGshhvEpM4R3i+Fc0rwyfHeTgVyDDONTZBDvvoU+vAzokrrO4D1bHEeufobxKTGIzy4DuqSuMHTPF7OBwDDeOP/Q5tNiH4iGdEl1MXivDLl7tl0FtcwwLtVo2IenYV1Sj2FbEhjGG9cbD+YRcsF4H74GeKn7DNYaV2yLtquglhnGp6T/JA2DuZaj7g95w71keFZ7PGlTgwzjLRj8QzSca5raCCF+AdBSDMaadwZwLcUw3gGGc807w5aklcTwreUwjHfQQn/EBnRJkrrH4K1JGcZnxGJ/7IZ0SZKaZ+hWUwzjM86QLklSPQzcaoNhfE4Ne0MxrEuSVhrDtrrIML5CjfqGZGiXJHWdIVuzzDCuJdX9Bme4lySBAVrqMYxrqtp+8/XLgCS1/14saQfDuFaUefgA8guF1L55eC+R1A2GcWnGGAIkSZofu7VdAUmSJGmlMoxLkiRJLTGMS5IkSS0xjEuSJEktMYxLkiRJLfFqKsuwahvsdXNOfbubnxxT36YkSRrNJNlglVerXfEM4zOgjS8Amh9+mZOG831WUlsM49KcM2RIktRdjhmXJElSZ0XEgRHx0Yi4JSK2RsSmiHhfROy7zHL2K6/bVMq5pZR7YNPbXopHxiVJktRJEbEWuAJ4PHARcD3wPOAM4MSIODoz7xqhnP1LOeuADcAFwHrgdOCkiDgyM/+3iW0P45FxSZIkddWHqcLwmzPzlMx8W2YeB5wDHAq8a8Ry3k0VxM/JzONLOadQBevHl+00te0lGcYlSZLUORFxMHACsAn40MDiM4H7gVMjYvWQclYDp5b1zxxY/MFS/ovL9mrd9igcprIMq7ZuZ58btzZS9r1r92ikXEmS1J5huWHV1u1TqslMOq5ML83MnToqMzdHxFepAvPzgS8tUc6RwJ6lnM0D5WyPiEuB1wHHAr2hKnVteyjD+OjW3P/AHfznlR9spvQrmylWkiR11/0P3AGwpuVqsOWhH3HFnf9Qe5nA+oi4aqHlmfnsIUUcWqbfWWT5d6kC8TqWDsSjlEMpp+5tD2UYH93PbN/+k4c3b7nl6rYr0iHry/T6VmvRLfbJzuyPXdknu7JPdmZ/7Gqe+2QNcF/Ldbh+Ow9x30N3NFH2mgleu0+Z3rvI8t78xzRQTl3bHsowPrprYaRvcStG75uufbKDfbIz+2NX9smu7JOd2R+7sk+alZm/0XYdxtS7q92kN9QYp5y6tu0JnJIkSeqk3tHnfRZZvvfAenWWU9e2hzKMS5IkqYtuKNN1iyw/pEwXG9c9STl1bXsow7gkSZK66LIyPSEidsqsEbEXcDTwY+BrQ8r5Wlnv6PK6/nJ2ozoRs397dW57KMO4JEmSOiczbwQupToJ9A0Di88GVgPnZ+b9vZkRsT4i1vevmJlbgE+U9c8aKOeNpfxL+u/AOc62x+UJnJIkSeqq36W6Jf0HIuJ44DrgCKprgn8H+KOB9a8r0xiY/3bgGOAPIuJw4OvAU4GXAreza+AeZ9tjicyJTwKVJEmSGhERBwHvBE4E9gduBf4FODsz7x5YNwEyczCMExH7Ud098xTgCcBdwOeBd2TmDybd9tjtM4xLkiRJ7XDMuCRJktQSw7gkSZLUEsO4JEmS1BLDuCRJktQSw7gkSZLUEsO4JEmS1BLD+ICIOCQi/jAiNkTE9yNiW0TcFhEXRcSxY5Z5VERcHBF3R8QDEfGtiPi9iFhVd/2bEBGPjIgzIuJjEfHN0icZEa8do6w15bWLPS5oog11qrM/+sqc6X2kp652DNlHJr71cJ0i4sCI+GhE3BIRWyNiU0S8LyL2XWY5+5XXbSrl3FLKPbCpujeljj6JiI1D9oNHNdmGukTEyyLiLyPiKxFxX6n7J8csq5Z9rW119Ulp/2L7xw+bqLvUBO/Auas/BV4B/DdwMXA3cChwMnByRJyRmR8YtbCIeCnwj8CDwGdLeb8MnAMcDby81to3YzXwvvLv24AfAgdNWObVVBfNH3TthOVOQ639MSf7SBPtuBk4b4H5C96YoQ0RsZbq7myPBy4CrgeeB5wBnBgRR2fmXSOUs38pZx2wAbgAWA+cDpwUEUf236a5y+rqkz5nLzL/oYkqOj1/DBwGbKHad9cvvfrCGujXNtXSJ8W97Hg/7rdlgjKl6cpMH30P4DTgmQvMfyGwDdgKPGHEsvamusXqVuA5ffMfRfWmmsAr227zCO3YHXhJr93AWaXurx2jrDXltee13a6O9Me87CO1tqOsv7Htdo1Qz0tKXd80MP+9Zf5HRiznr8v67x2Y/+Yy/wttt7WFPtlYfUS136YJ++NY4BCqW3MfU/rgk231axceNfbJJmBT2+3x4WPSh8NUBmTmeZn5jQXmf5nqw2F34KgRi3sZ8Djggsy8sq+sB6mODAD8zkQVnoLM3JaZn8/MW9uuSxfU3B9zsY8wP+0YWUQcDJxAFQg+NLD4TOB+4NSIWD2knNXAqWX9MwcWf7CU/+KyvU6rq0/mSWZelpnfzcyxb3c9b/1aR59I88RhKsvzkzId9efR48r0Cwssuxx4ADgqIvbIzK2TVm7G/GxEvB7YH7gL+PfM/FbLdWrDvOwjTbTjMRHxm8ABVD9FX5WZXRov3mvzpZm5vX9BZm6OiK9SBajnA19aopwjgT1LOZsHytkeEZcCr6M6mtj1oSp19clPRcQrgKdQ/TJ5HbCh438LTai9X+fIHhHxKuBJVF9KvgVcnpkPt1staXSG8RFFxJOB46lCxeUjvuzQMv3O4ILMfCgibgKeDhxM9SGzkryoPH4qIjYCr8nM77VSo3bMyz7SRDsOA/6uf0ZEXA2cmpnXTFDXuiza5uK7VAFpHUsHpFHKoZTTdXX1Sb/Bk7pvj4g3ZOaFY9RvVjXRr/PiAOATA/NuiojTyy/aUuc5TGUEEbEH8ClgD+CszPzRiC/dp0zvXWR5b/5jJqjerHmA6iTZZwP7lscLgcuoxg5+aVZ+aq3JvOwjdbfjvVQnfT4O2At4LnAhVUDfEBFPHLOedaqrzfOyD0C9bbmI6gTgA6l+OVgP/Hl57Wcj4iUT1HPWzNM+UqePUR0kO4DqxPpfoDr/Yg3w+Yg4rL2qSaObyzA+5HJHCz0WvaRSuSTbJ6iCwWeB99RZ1TJtfNxcnX0yicy8PTPfkZn/lZn3lMflVEd1/gP4OWDsSwSOqiv9MUpVy3TW95FltSMz35KZV2TmnZm5JTOvzMyXU12t5bHAW5fZvDbU9X83tX1gCkZuS2aek5mfy8z/y8wHM/OGzHw78Baqz653N1nRGTNP+8jIMvPszNyQmbdl5gOZeW1m/jbVl/k9qU6ulzpvXoep3Eh1ebVR3bLQzBLEP0l1Sba/B161zBNOekcr9llk+d4D6zWplj5pShnKcC5wBPAC4P0Nb7Ir/TEv+8i02vER4Neo9pG21dXmLu0Dk5pGW86lulzm4RGx1+A4+zk1T/vINHyE6ktbF94npKHmMoxn5vGTlhERjwA+TRXEPw28eowTQm4AnkM1ju+qBcp/CtXJoI2flFVHn0zBHWXa+DCVDvXHvOwj02rH1PaREdxQpouN5T6kTBcb51t3OV3QeFsy88GI2Ew1xG01sBLC+DztI9Nwe5l24X1CGmouh6lMKiJ2pxqf+nLgfKoTxsY5M3tDmZ64wLIXAI8GrliBVwZYzPPLtOtXjKjTvOwj02pHl/aRy8r0hIjY6b00IvaiGtr2Y2DYFWC+VtY7uryuv5zdqIZw9W+vy+rqk0VFxKFUQXwzcOe45cyYxvt1zhxZpl14n5CGMowPKCdr/jPwUqorOZw+eCmpBV6zT0Ssj4gnDCy6kOrD4pUR8Zy+9R8F/Fl5+le1Vb5DFuuTiDiifNkZXP844PfL07bGZzdmBewjy25HRDy69MmTBuY/a6GTeCPiGcC7ytPW95HMvBG4lOpksTcMLD6b6qjc+Zl5f29mae9OdxvMzC1U56WsZtcxrm8s5V+SM3AHzrr6JCIOXugk3Yh4LNVJe1Bd035W7sI5koh4ZOmPtf3zx+nXebFYn0TE0yNivwXWfzLV9fmhA+8T0ijCa+7vLCI+RnUXzjuBD7PwCTEbM3Nj32tOo/qA+HhmnjZQ3ilUQeVBqkt03Q2cTHWpqguBX5+FGx9ExNvYccviw6muanEFOy679m+ZeW7f+qexQJ+Uyxc+neoGSr3bmj+DHdfR/ZPM7IW3zqqrP8qyedlHltWOiDiG6ojflzPzmL755wG/SnW0/ftUd/VcT3XUfRXwt8Dru9Ansestyq+jOu/hWKohA0dl3y3KI6K6vWhmDJSzfylnHVW7vw48leqgwO2lnBubbk8d6uiT8vdyLvBlqnMZ7qa6jvQvUY2bvhJ4UWbe03yLJlP+Lk4pTw8AXkx1xPYrZd6dmfnWsu4a4Cbg5sxcM1DOsvq1y+rok4g4C3gb1XvITVS/lKwFTqK68+/FwK9k5rZGGyPVIRu+xeesPSi3YB7yOGvgNaexxC3eqX5CvBj4EdVPiddQHQVe1XZ7a+yX8wbWX7BPgN8CPkd1J7ktVEHre1RXqvnFtts57f6Yp31kue1gx22wNw7MPwX4J+B/gPuobvZyK/CvwMltt3GBdhxE9UXr1lLXm6lOQN5vgXWTRW7xDuxXXndzX5s/ChzYdhun3SdUl6g7r+w/d1HdcO1uqrD2JmD3ttu4jL44a8h7xaa+ddcMzhu3X7v8qKNPqC6J+xngeuCeso/cAXwReDXlYKMPH7Pw8Mi4JEmS1BLHjEuSJEktMYxLkiRJLTGMS5IkSS0xjEuSJEktMYxLkiRJLTGMS5IkSS0xjEuSJEktMYxLkiRJLTGMS5IkSS0xjEuSJEktMYxLkiRJLTGMS5IkSS0xjEuSJEktMYxLkiRJLTGMS9KEIuLciMiI+OIS6+wVEdvLeidNs36SpO4yjEvS5J5Rpt9cYp3DgSj/vrrZ6kiSZoVhXJImEBG7AT9fnn5jiVUPL9M7M/MHzdZKkjQrDOOSNJl1wJ7l38OOjINHxSVJfQzjkjSZ3hCVHwM3LLHeM8t0qcAuSVphDOOSNJnDyvSazHx4oRUi4hHA08pTw7gk6acM45I0mV4YXypkPw3YY4T1JEkrjGFckibTG6ay1MmbvSEqDwLXN1sdSdIsMYxL0pgiYl/goPJ0qSPevTD+7cx8qNlaSZJmiWFcksZ3WN+/r1livWPK1CEqkqSdGMYlaXy9ISoPUw1B2UVEPJfRxpVLklYgw7gkja8XslcBRwwujIi9gb/pm+U1xiVJOzGMS9L4+oepfCoiXhERayPikIg4HbiKHTf7AVgbEUdNt4qSpC6LzGy7DpI0cyJiFbCZ6u6b5wKvXWC1+4D3AO/sm3d+Zr6m+RpKkmaBR8YlaTzrqII4wNnAqVQncW4FbgM+DTwL+AvgEqo7dF4JfHzqNZUkdZZHxiVpDBHxSuAzwD2ZuW/b9ZEkzSaPjEvSeHpXUvl2q7WQJM00w7gkjad38ua1rdZCkjTTDOOSNJ5eGPfIuCRpbIZxSVqmiNgPeGJ56pFxSdLYDOOStHz91xc3jEuSxubVVCRJkqSWeGRckiRJaolhXJIkSWqJYVySJElqiWFckiRJaolhXJIkSWqJYVySJElqiWFckiRJaolhXJIkSWqJYVySJElqiWFckiRJaolhXJIkSWqJYVySJElqiWFckiRJaolhXJIkSWqJYVySJElqyf8D2HGlAdCIA/wAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 282,
       "width": 369
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def norm_gamma(x, mu, lam, a, b):\n",
    "    mux, taux = x\n",
    "    \n",
    "    pmu = norm.pdf(mux, loc=mu, scale=1 / (taux * lam))\n",
    "    ptau = gamma.pdf(taux, a=a, scale=b)\n",
    "    \n",
    "    return pmu * ptau\n",
    "\n",
    "a, b = 1, 1\n",
    "mu, tau = 0, 1\n",
    "\n",
    "X = np.mgrid[-2:2:0.1, 0.1:5:0.1]\n",
    "Z = np.apply_along_axis(norm_gamma, 0, X,\n",
    "                        mu=mu, lam=tau, a=a, b=b)\n",
    "\n",
    "plt.title(\"Normal-Gamma\", fontsize=14)\n",
    "plt.contourf(*X, Z)\n",
    "plt.xlabel(r\"$\\mu$\", fontsize=13)\n",
    "plt.ylabel(r\"$\\tau$\", fontsize=13)\n",
    "plt.colorbar();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Posterior Distribution with Improper Priors\n",
    "\n",
    "In order to simplify the estimation of posterior parameters, we consider improper priors such that $\\mu_0 = \\lambda_0 = a_0 = b_0 = 0$. \n",
    "\n",
    "Note that $\\text{Gam}(a, b)$, in scipy, is defined as:\n",
    "\n",
    "$$\n",
    "    f(y, a) / b = \\frac{x^{a-1}\\exp(-x/b)}{b^a\\Gamma(a)}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 159,
   "metadata": {},
   "outputs": [],
   "source": [
    "E_tau = (N - 1) / ((D - D.mean()) ** 2).sum()\n",
    "E_mu = (D ** 2).sum() - N * D.mean() ** 2 + N / E_tau\n",
    "\n",
    "def q_mu(mu, D, E_tau):\n",
    "    N = len(D)\n",
    "    mean = D.mean()\n",
    "    precision = np.sqrt(1 / N * E_tau)\n",
    "    return norm.pdf(mu, loc=mean, scale=precision)\n",
    "\n",
    "def q_tau(tau, D, E_mu):\n",
    "    N = len(D)\n",
    "    return gamma.pdf(tau, a=N/2, scale=2 / E_mu)\n",
    "\n",
    "\n",
    "def q(x, D, E_mu, E_tau):\n",
    "    mu, tau = x\n",
    "    qmu = q_mu(mu, D, E_tau)\n",
    "    qtau = q_tau(tau, D, E_mu)\n",
    "    \n",
    "    qx = qmu * qtau\n",
    "    return qx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 158,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwoAAAI/CAYAAADX+C6LAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAWJQAAFiUBSVIk8AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd9wsZX3//9eHJu1QFRNBpEgxivoVCyUqJUEwsWEJFhSJRqOxJCa/GE0imJiuIH41sQSPxlhiw29sISoQBYlKbBiayEEjJhBAOMChnuv3xzXrvWdn+87szs79ej4e+9i9d9q1187OPe+5rpmJlBKSJEmS1G2zRRdAkiRJUvMYFCRJkiSVGBQkSZIklRgUJEmSJJUYFCRJkiSVGBQkSZIklRgUJEmSJJUYFCRJkiSVGBQkSZIklRgUJEmSJJUYFCRJkiSVGBQkSZIklRgUJEmSJJUYFCRJkiSVGBQkSZIklRgUJEmSJhARJ0XEuYsux7KzHpvPoCBJkjRCRGw5y3Bl1uNyMShIkiQNEREPBa6IiKf0GbZnRPw/4C/mX7LlYj0uH4OCJEnScJcAbwXeGxEfB+4HbB4RrwYuBm4Fzlhg+ZaF9bhkDAqSNIWI2CsiUvFYu+jy1K3oS9z5vCf1Gb7U9TGq/BFxRNfwU+Zfwmos+/e0KCmlu1JKpwEPBH4MnAL8IvBrwDEppWenlK6eZJ4R8bTie7g9InavvNCbLuvEYlk/jYjd6lzWMHXUI8ynLptSh/NmUNBIXf9Uxn2cvugyD1L8kzyleByx6PKozO9I/bhe9Ge9zF0CNnb9vbHn77FExNbAW4o/35VS+nEFZRvmg8DlwI7An9e8rHFUUo8w17psWh3OhUFBq81ewBuKxxELLYkG2Qu/I5XthetFP3thvdQuIraIiFcA3wf2JB8J/wrwMeALEfGBiNhjglm+jPzd3c4c+uSnlO4B/rT486SIeFDdy+ynhnqEOdVlU+pw3rZYdAG0dJ42xjhX1l4KacFSSuuAWHQ5mmLZ62PZyz+u1fI5a/Bg4HeBk1NKnyq6392TUjqt6Gv/duB3isdQEbEN8Nriz7UppWtqKnOvDwJ/AjyAHCxPmNNyu1VWj7CQumxCHc6VQUETSSmdtegySJI0Tymlb0fEA1NKd/UZ9kPgSRNc1vP5wH2K1++vqoyjpJTuiYh/BF4HPCMi9izKPjcV1yPMuS6bUIfzZtcjSZKkEfrt3E4yvMtvFs9XppS+OlupJvaB4nlz4DfmvGyg0nqExdTlwutwrlJKPnwMfZBPOkp5dZl6HjsAzwX+Hvgm8FPgLuAG4D+ANwP7TjC/bcn9Ej8N/AjYUDx+AHyC/OPdoWv8I7o/x7DHgOU9iHxJt4uBm4plXQ38E/C0McrbvfxTivcOAE4nXy7u5mLYSRPWa7/5HgS8i9wFbANwHfAF4NkTzHfWz7s5cCLwz8X3c3sxjx8V3/e7gOOBbSv8jrYHXg38K3ANcEexfn0deCNwnyq/I3Kf2M74a+us07rWn55l/Grxff138X2tA/4ROLQYflJXGUrLGac+5rleVP19DpjfVL81cr/szryOmGbcGepl5PdU4XagX53tSd7eX0q+FOZPgQvI2/Mtpl1/Byx/V+CPgYuK8t8MfAf4Q2D7Ypyri/J9t8plDynTQV118sYxxr9P1/hnjhh3O+CeYtwPDRnvm8U4VwHR9rokn3x8D2P+Xvo8dq2iDpf1sfAC+Gj+Y9A/nAmm34q8QzDqx3gX8PIx5ncseWdm1Pze2zXNEWOM3/czAqcCd4+Y7tx+G5MByz+F3Fx6W5/5nDRh3fbO98QRdf1pYOsR85zp8wL3Br42Zn0/taLv6Djgf0ZMdzPw5Kq+Iybb4Zq1TmtZf4p5bw6sHVKue4D/jxmDwrzXi6q/zz7zm/q3xhIEBerZ7h0L3DhkfmcD95p0HR6w7GeNWNbXgZ8jX2knAe+vYrljlOuPuspw5BjjH9M1/itGjHtY17i/P2S8t3SN9/C212Wx3o31W+nzuKyKOlzmh+coaB42A+5FPsr7r+SjEP9D3qjcn7xxexL5nJn/GxHXpJQ+2W9GEfEs8slEmxdvfQf4OPmIXvf8nsCmJ+xdTD4R+yHkE5EAPgJ8eFjBI+LPWTlR6p5i/C+Rj6wdBJwM3Bd4PPCliDgkpbRhaG3A4cDri/n9PXA+eYfjAHIAmtajyP0mAc4E/q1YxqOAXycfbfoVcrPpM/rNoKLP++5imZCvbPEh8iXlNpBblg4AHgc8pme6ab+jpxfjbV6U+dPAF8l1uQY4knyN7jXAJyPil1NKXxo2Tyr8jmpYh6pef84AXlC8vhN4H/kqJBuBR5PXnb8EZj0/aa7rRY+q62zm31oFqqiXvmra7j0c+D3ydvmdwFfJrX6PBF5KrrNfJn9Pfzxj+V8AvLdY1l3AJ8nbhJvIJ9O+rFjuu1n5P/HNWZY5gV8unjcC3xhj/P/T9XpUGccd98Ku108AvjVoxJbU5eX0vxDL48mt0JBPov5Cn3EGbR/GrsOlt+ik4qP5DwYcmZpg+s3JiX5g8xzwUOAnxXKuBDbrM87ewC2sHOV81aB5AjvT50gdfZrCh5TpUFaOkNwCPK7POLuQj6Z05vnXA+bVvdxUfNZfqOC76Z3vzcAhfcbbj3xzm854T6/j8wK7sdLE+3VguyFlfwDwgBm/o/uT/2El8gb9UQPGexS5i0Mid3PZctbviPG62lSyDtW4/jy2q3w3Agf3GecAVn6bncdJk9THvNeLOr7PPvOb5bd2StfwI0Z8jqHjTlEvoz5nndu9q4H9+oz3aPJOaCJ3F5y6VYG8s7yhmNc6+hztLerg5p6yDf0eqniQ/xfeWizv4jGn+XAx/kZgzYhx/77r8+w2ZLwHdI33idVSl33m8addZX7shNOOVYdteHgysyYyxs3W1vZOk1K6J6X0+VT8uvpJKX2HlSN0+5BbBXr9AfmoE8BfppTeOmieKaUbU0rnTvLZ+ugc/QL4vZTSv/VZzg3kI4a3FW/9ZkTsNMa8X5JS+s8Zy9fP76WULux9M6V0BflIZ8fv9puW2T/vPqxcJOGDKaVbBxU0pXR1muIOnH3KvEPx+pkppa8PWNbXWbnc3h7AM8eYdxXfUV3rUFXrz2u6yvfqlNJFfcp3GfCiGZcz7/Win6p/c7P81pqszu3e84r66Z3f18itIZAP8jx64lIDEbE5+byarYH15Dv9lo70pnyJ2Lf2vD2PI8L7kM+vA7hszGk6rQRXpJTWjznuNSmlaweNVPy+Oi1AD+03Tkvrslf3Z//OJBOOU4dtYVBQk1zQ9XqT7gfFRuvXij/XU/NdESPiXsATiz+vJx+p6avYYHyo+HM7cp/SYa4mn8xZtRvJTcR9pZQ+D3R2lA6JiJ/rDKvw897W9frBY5d8ChER5BPkAb6WUvryiEk+Qu5zDXP4jmpchypZf4ryHVf8eS0rV/LoV77PkE8Antbc1osBqv7NTf1ba7Kat3vfHPEb7e4O+Asj5jXIc8knYENuXbl8yLhf6Xp9VUrpp1MucxIP6Hp9w6iRI2J74IHFn0O78xSXFO38tsbZUb+xeL5/sS3t1aq6HOBhxfMPU0o3TTH9qDpsBc9R0KRG3XBt4PWEI2Ivcl/oI4ADgZ3IRyv66b0z40NZOXJ8zhhHVmb1MPJ5FQDnppTuHDH+2awcRXwM+aogg3xlWOvKDL48Rjm/xMo/4UexsvNU1ef9HvlclPsBv15sPN9N3pHfONanGN+DyV0gAG6IiKeOMc0t5PVu1B01q/iO6lqHqlp/Hka+0ECnfPeMGP+LjK63Qea5XvRT9W9ult9ak9W53Su1vvT4cdfrnUeMO8hriuebgL8bMe5/db2eV5/6Xbpej7Nz+zBWDuiOKuNDWPk9j/N5rif/HrciB71beoa3rS43UbSA7Vn8+e0pyzCqDlvBoKCJpClvuBYRrybfWv1eo8Yt7NDzd3dwmOXI5rh+vuv1sCMp/cb5+YFjZT8eMXxa359wnPt1va7k86Z8M5qXkE8w34p80uPJwE8j4qvkI0//0q+LyxT26np9bPEY16gdkSq+o7rWoarWn+7vf9J1ZyJzXi/6qfo3N8tvrcnq3O7974jhd3S9HnQAaaCI2JOVLiCfSCndNmx8VrqtwPxORO3+/zfOwa5pT2Qe5/Pc3PV6G7p2cltal72m7nbUZWAdtoldj1S7iHgucBorP+wvk08iehG5O9HTisdLuibbnE11B4d5/BjXdL0e2J+6S3eZ1gwcKxt1dZBpjdqYw6afZfuu15V93pTSp8l9jM8in6AI+Sj+ccCbgG9ExHcjYpId+352nGHarUYMr+I7qmsdqmr96f7+J113JjbH9aKfqn9zs/zWmqzO7V7dLUe/3PW6dF5FH91dV+Z1FLw7DPUeDOunjisedXRvP3t/H22sy14P63o9bVAYVoetYYuC5uGNxfPd5OvYf67fSBExrO9yd3Kfxz/d7iMU2w0ca0V3meruFjXItqNH2eSzdP+Tr/TzppS+DTwtItaQL015GPnSl4cBW5KbyT8bESemlP5xjOX1013+U1JKp045n7o0fR3qrr9J152pzGm9mIdZfmuTmucBvaavs8Ps2fX6yjHG/8Wu1/Paue3uIrPLwLFWdHb+/yuldN2IcR9RPK8n33h0lM7y76QcCttYl72qCArD6rA1bFFQrSJiH/LVCQDOGhQSCg8YMqy7D+S0/aQn8ZOu1/uNMX73ONdUXJZxPXD0KJuM013OWj5vSml9ccWrP04pHUHunnBaMTiAtxQnqk+juzvJIk6QHaXp61D3MiZdd2ZS83oxD7P81mDTo6GjWrfuPVaJqtH0dXaY+3a9HtriU5z4e0Lx53Uppbq6g/Za1/V66M5tz8nJQ/vQR8Q2wMHFn98a83yczvJ/2Gf8VtXlAAcVz7cDpStxjWlYHbaGQUF1697gjDoy8YQhw77DSqvCkcURyWl0N38Pu0rBt1n5Z35EsTEcpvuKH1+bpmAV+MWIGLXTcWTX6+5Lic7l86aUrk8p/Q4rN8fZjfIOybjf0TdZWSeOiYiZj3hXrOnr0LfJR8IAHj/GjvlRdRWk4vViHmb5rUG+p0fHqPMXem9A16vKemn6OjtM98n4u48Y9+Ws/G+a1xFwgKtYOfJ8wIhxH8xKiBx1U8CnsdK1d+Q5AsWFRTrngfQ7mt62uuxn/+L58jEu5FAyRh22hkFBdes+GrHvoJEiYg/ghYOGFz/kzqX41pDvqTCN7i4AA3csU0p3AJ8p/rw3cNKgcSPi/sCziz9vJV8JZBF2YeUOuyURcQwrR6i+mlL62T+fBXzedV2ve7tAjvsd3UO+zjfkvqKvGzTuIjR9HSrK99niz/sCzxk0bkQcx/SXrJzEuq7XU60XczL1b63QfT+HgQEsIg5npUvJIJXVS9PX2RG6773xi4NGioiHkc+R65jbzm2xzeqcsH9gRAzrW999zsHAnfWiNeENXW+N83m6w+e/9xnetrrcRBHyO/f9mLYFZFQdtoZBQXW7hJXU/5SIKN1IJyLuC3yK0SfD/SUr/xR/PyJeNejaxRGxU0Q8vs+gq7pej/oH/NesHK17c/FPu3c5OwMfY+Uf9N/O6RrSg/xNRDyq982I2Bc4s+utN/eZdubPGxFPKL6XgScaR8QDWTlZ7hbKLU2TfEd/xsrR2T+IiN+NiIHbtYi4T0T8YUTM6wY5TV+HuteDt0bEw3tHiIj9GHI9/XEsYL2Yh1l+axeyst6eEBEH945QzGfgvS26VF0vTV9nB+m+D8NvRkSpK2vxWT7PpoFqnkfBAf61eN4MeOSQ8bqDwpHFTvkmit/Tx1k5Og7j3Xyseyf3X/oMb1td9tPpKjTtOY+j6rA1PJlZtUop3RkR7yTfFXdL4N8i4kxyU/xd5H9sLySn+/cDzx8yr6si4tfJLQubAacDJ0fEx8iXIkzkIy+Hkq+k8lHgvJ553BgR3yRvhI+MiL8jXx9+fdc4ny+eL4yIvyS3XqwBzouID5E3ohvIJ16+iJVm1+8AfzxFNVXls+QdrfMj4n3kq0vdQ76G+6+zskH8RErp470TV/R5f578vfxVRJxDPtLyA3LL0r2LsjyLlX8up6eUNrlaxITf0X9FxAnA/yM30/818BsR8XFySL2NfEWM/YBDgMeSr6h17rCKrErT16GU0lci4h3Ay8iXjL2wWHe+Qt5ZfDR53dmOfLWice5V0c9c14s5mPW3dkdEvA34I/J28dzi83yD3IXkUPK2MMjr9pMHFaTqemn6OjtISukbEXE++ST57YGvF+v2f5JbgJ4APIlcL+eS7+cDA3Zui1a0zwLvSSm9uM/wN5P/rz2j33c8xCeBPyleH8GmO+XdOkHhDvI6cUFErAUuJv82DyKfG7Ar+TN1DrT9RkQ8GfijlNJd9NfpFndVcYGBTVRZl0WXxpsYr7VrbUppYM+CPsaty00U+yXrgL3J3QhPJ9/w9U7gc0XL2ihD67BVUko+fAx9kHfAU15dppp+a/IPOA15/B35pOfO32uHzO9XgetGzC8BZw6Y/jjyFZj6Ttdn/DcOG794nAvsOqTMR3SNe0qF380m8wWeRz45a1A5PwNsPWKeU39e8s7NqO8lkf/RnQ5sVtF3dAj5CPQ4y14PHDTrd0S+j8PI9bWKdaiu9aeY9+bA+4aU6x7g98jdUDrvnTRJfSxivaj6++ydHzP+1sg7f2cPmf4m8rbulK73jqigXsZab2ddZyf9Dib9vgbMY1/yhS8GlfcK8pHn84q/bxmyrr2+GOdlA4Z3/qftO0U5v1lM+/0Bw4N8/lUCvkBuWRr0mU4B3tPz3iVDlv2grvH+pO66BH4OWNvzuJP8W+99/5iq63LIdL/d5zPdOmh9mKYO2/JYeAF8NP/R/UOaYR5bAL9JTu03k//BrgM+0tk4MNmO15rih/4F8oled5KPTn6f3JJwErDdkOkfQ+7f3jmqOfQzkvtnn0G+w2yn/D8iN78fP8bnP6J7w17hd1OaL/lI07vJO84byHeP/CLwnAnmO9XnJf+De2jx3ZxV/DO5hbzD8dNio/424P+MUYZJv6MtgOcW69QPyIHgruLzfx14F8VR6yq+o0nW11nXobrWn55l/CrwaeDaomxXAx8EDi2Gn9RVhpMmqY9FrBdVf5/95seMv7VinX0pcD45GGwo6uatwN7FOKd0LfeICupl7PV2lnV20u+gqnWcfCL86UU93FGsz18GXgFsW4zz02I55w+Zz8eKcQ4dMPyGok5iijK+uOuzHtZn+H5dw/+mWE9eB1xafAfXkn+rnf+d+5L7699KbpU6Yciy31TM925gz3nUZc88dy7Gv2La73iSuhwx7Uu66i0BF4453dh12IZHFB9a0hKKiCOAc4o/T00pnbK40khSsxWX7O6c//L2lNJvDRjvB+RLdu+QUrq1Z9je5J3nr6SUHjtFGbYmB/HdgHellF7SM/xZ5AMeAM9LFd1TpOgG9H1yWPxISumE4VOMnN9YddkzzePJLVGfTCkdP8vyi/kNrcuqVV2Hy8CTmSVJ0mrRfcL+oPMTdiT3X7+8NyQUOucPjLwUaT8ppduBvyj+fH5E9F4id9K7LI/r2eQd3I1AFTeoHFmXfXROyq7kkqJj1GXVqq7DxjMoSJKk1WKcnfBHjBg+U1Ao/C35ilVbU77cd2f+GxjvKkYjFUfC/7D4c21K6ZIKZjtNoOlcce67FSy/Y1hdVqamOmw8g4IkSVotOkfB7yKfe9FPZwd4VJCYOigUR8J/p/jzxRHRfa+EzvK/m6a4GdgAzybfmOwmqtuZHqcue1XaogAj67JKddRh4xkUJEnSatHZCf/PNPgymJ0gUDrqHRH3Ah5HPpH14lkKklI6K6UUKaWtU0o/LuZ/P3J/e6iw21FK6QPFsnZKKV1b0WzHqcufKe579AvklpLe+6TMpF9dVq2mOmw8g4IkSWq9iNiVlbscD2sN6OwA39Zn2LPJ9xa4dJyd4yl0d+eZpWtTrSaoy267AdsCV6aUNo4aWc3gDdekJZZSOpd86UlJ0hAppesZsb2MiG3J3UsAnhcRX07F5SEj4mjy5UIBLq+pjJ8ZVcYmGKcu++jsc24ZEdGpVzWbLQqSJEnZQ8k3IfwB+Rr9l0XEpyPiYvJ9ezonF/9SRHxiQWVcVj8h3/foAODfI+J9EXHvBZdJIxgUJEmSss75Cf9IDgpbAr9EvsHWycDTyf3r7wb+YxEFXFZFd6MTyF2VHg48g3zjOjWYN1yTJEkCIuLdwIuAp6WUzlp0eaRFs0VBkiQp67QoXLTQUkgNYYuCJEla9SJiS2A9sD6ldJ9Fl0dqAlsUJEmS8jX+74XnHkg/Y4uCJEmSpBJbFCRJkiSVeMO1BYiIq4AdgHULLookSZLabS/g5pTS3pNOaFBYjB02Y4tdtt9y510WXRBJkiS11y133chG7p5qWoPCYqzbfsuddzlstxMWXQ5JkiS12AXXfpib77pu3TTTeo6CJEmSpBKDgiRJkqQSg4IkSZKkEoOCJEmSpBKDgiRJkqQSg4IkSZKkEoOCJEmSpBKDgiRJkqQSg4IkSZKkEoOCJEmSpBKDgiRJkqQSg4IkSZKkEoOCJEmSpBKDgiRJkqQSg4IkSZKkEoOCJEmSpBKDgiRJkqQSg4IkSZKkEoOCJEmSpBKDgiRJkqQSg4IkSZKkEoOCJEmSpBKDgiRJkqQSg4IkSZKkEoOCJEmSpBKDgiRJkqQSg4IkSZKkEoOCJEmSpBKDgiRJkqQSg4IkSZKkEoOCJEmSpBKDgiRJkqQSg4IkSZKkEoOCJEmSpBKDgiRJkqQSg4IkSZKkEoOCJEmSpBKDgiRJkqQSg4IkSZKkEoOCJEmSpBKDgiRJkqQSg4IkSZKkEoOCJEmSpBKDgiRJkqQSg4IkSZKkEoOCJEmSpBKDgiRJkqQSg4IkSZKkEoOCJEmSpBKDgiRJkqQSg4IkSZKkEoOCJEmSpBKDgiRJkqQSg4IkSZKkEoOCJEmSpBKDgiRJkqQSg4IkSZKkEoOCJEmSpBKDgiRJkqQSg4IkSZKkEoOCJEmSpBKDgiRJkqQSg4IkSZKkEoOCJEmSpBKDgiRJkqSS1gWFiNgjIs6MiGsi4o6IWBcRp0fEzhPO5xcj4lPF9LdHxA8j4rMRcWxdZZckSZKaolVBISL2BS4CXgh8DTgN+AHwKuCrEbHrmPP5TeDLwNHF82nAecDjgc9FxOurL70kSZLUHFssugAVewewG/DKlNLbOm9GxFuA3wbeBLx02AwiYkvgz4HbgYNTSpd1Dfsz4JvA6yPib1JKd1T/ESRJkqTFa02LQkTsAxwDrAPe3jP4DcCtwIkRsd2IWe0C7Ahc3h0SAFJKlwCXA9sA21dQbEmSJKmRWhMUgKOK57NTShu7B6SU1gPnA9sCh4yYz7XAdcD+EbFf94CI2B/YD/hWSun6SkotSZIkNVCbuh4dUDxfPmD4FeQWh/2BLw6aSUopRcTLgQ8AF0XEJ4FrgN2BpwHfA04Yp0ARcdGAQQeOM70kSZK0KG0KCjsWzzcNGN55f6dRM0opfTQirgE+BDy/a9D/AO8lnyAtSZIktVabuh6NEsVzGjlixPOAL5CvePQgcpelB5FbIv4v8OFxFphSOrjfA7h0mg8gSZIkzUubgkKnxWDHAcN36Bmvr+I8hDPJXYxOTCldmlLakFK6FDiRfPnVZ0bEEbMXWZIkSWqmNgWFzhWK9h8wvHNi8qBzGDqOAbYEzutzUvRG4N+KPw+eppCSJEnSMmhTUDineD4mIjb5XBGxBjgc2ABcOGI+9yqe7zNgeOf9O6cppCRJkrQMWhMUUkpXAmcDewEv7xl8KrAd8P6U0q2dNyPiwIjovQLRl4vnZ0TEQ7sHRMTDgWeQz3P4UnWllyRJkpqlTVc9AngZcAFwRkQcDVwCPAY4ktzl6PU9419SPHdOdCal9LWIeC/wQuDrxeVRryYHkKcCWwGnp5S+V+PnkCRJkhaqVUEhpXRlRDwSeCNwLPBE4CfAGcCpKaUbxpzVr5PPRTgJeAKwBrgZ+Arw7pTSWFc9kiRJkpZVq4ICQErpR+TWgHHGjQHvJ2Bt8ZAkSZJWndacoyBJkiSpOgYFSZIkSSUGBUmSJEklBgVJkiRJJQYFSZIkSSUGBUmSJEklBgVJkiRJJQYFSZIkSSUGBUmSJEklBgVJkiRJJQYFSZIkSSUGBUmSJEklBgVJkiRJJQYFSZIkSSUGBUmSJEklBgVJkiRJJQYFSZIkSSUGBUmSJEklBgVJkiRJJQYFSZIkSSUGBUmSJEklBgVJkiRJJQYFSZIkSSUGBUmSJEklBgVJkiRJJQYFSZIkSSUGBUmSJEklBgVJkiRJJQYFSZIkSSUGBUmSJEklBgVJkiRJJQYFSZIkSSUGBUmSJEklBgVJkiRJJQYFSZIkSSUGBUmSJEklBgVJkiRJJQYFSZIkSSUGBUmSJEklBgVJkiRJJQYFSZIkSSUGBUmSJEklBgVJkiRJJQYFSZIkSSUGBUmSJEklBgVJkiRJJQYFSZIkSSUGBUmSJEklBgVJkiRJJQYFSZIkSSUGBUmSJEklBgVJkiRJJQYFSZIkSSUGBUmSJEklBgVJkiRJJQYFSZIkSSUGBUmSJEklBgVJkiRJJQYFSZIkSSUGBUmSJEklBgVJkiRJJQYFSZIkSSUGBUmSJEklBgVJkiRJJQYFSZIkSSUGBUmSJEklBgVJkiRJJQYFSZIkSSUGBUmSJEklBgVJkiRJJQYFSZIkSSUGBUmSJEklBgVJkiRJJQYFSZIkSSUGBUmSJEklBgVJkiRJJa0LChGxR0ScGRHXRMQdEbEuIk6PiJ2nmNdBEfH+iPhRMa9rI+K8iHh+HWWXJEmSmmKLRRegShGxL3ABsBvwKeBS4NHAq4BjI+LwlNL1Y87rJOA9wG3Ap4F1wE7AQ4AnAu+vuPiSJElSY7QqKADvIIeEV6aU3tZ5MyLeAvw28CbgpaNmEhGHkEPCxcCxKaX/7hm+ZZWFliRJkpqmNV2PImIf4Bjykf+39wx+A3ArcGJEbDfG7P4K2Bx4Xm9IAEgp3TVbaSVJkqRma1OLwlHF89kppY3dA1JK6yPifHKQOKnWAZsAACAASURBVAT44qCZRMQewGOBbwDfi4gjgYOBBHwLOKd3/pIkSVLbtCkoHFA8Xz5g+BXkoLA/Q4IC8Kiu8b8EHNEz/LsRcXxK6ftTllOSJElqvDYFhR2L55sGDO+8v9OI+exWPD8L+F/geHKwuA+5C9OJwGci4qCU0p3DZhQRFw0YdOCIMkiSJEkL1ZpzFMYQxXMaMd7mXc8vSil9MqV0c0rpSuAF5C5J+wNPr6eYkiRJ0uK1qUWh02Kw44DhO/SMN8iNxfMdwGe7B6SUUkR8Cngk+bKrHxo2o5TSwf3eL1oaHjGiHJIkSdLCtKlF4bLief8Bw/crngedw9A7n/UDTlruBIltJiibJEmStFTaFBTOKZ6PiYhNPldErAEOBzYAF46Yz3fI5ybcOyLu22f4Q4rnddMXVZIkSWq21gSF4hyCs4G9gJf3DD4V2A54f0rp1s6bEXFgRGxyYnFK6W7gncWff9UdOiLiIOAk4G7gYxV/BEmSJKkx2nSOAsDLgAuAMyLiaOAS4DHAkeQuR6/vGf+S4jl63v8z4Gjg+cBBEXEu+apHTwe2Bl7j5VElSZLUZq1pUYCftSo8ElhLDgivAfYFzgAOTSldP+Z8biMHhVOBbcktFE8mh5AnppTeUnnhJUmSpAZpW4sCKaUfAS8cc9zeloTuYbcBpxQPSZIkaVVpVYuCJEmSpGoYFCRJkiSVGBQkSZIklRgUJEmSJJUYFCRJkiSVGBQkSZIklRgUJEmSJJUYFCRJkiSVGBQkSZIklRgUJEmSJJUYFCRJkiSVGBQkSZIklRgUJEmSJJUYFCRJkiSVGBQkSZIklRgUJEmSJJUYFCRJkiSVGBQkSZIklRgUJEmSJJUYFCRJkiSVGBQkSZIklRgUJEmSJJUYFCRJkiSVGBQkSZIklRgUJEmSJJUYFCRJkiSVGBQkSZIklRgUJEmSJJUYFCRJkiSVGBQkSZIklRgUJEmSJJUYFCRJkiSVGBQkSZIklRgUJEmSJJVssegCSJKkam14yO5zXd42F/94rsuTNB8GBUmSltS8A8Ego8phkJCWk0FBkqSGa0ogmFZv+Q0O0nIwKEiS1DDLHgxG6f58hgapuSoLChHxAeCClNI7qpqnJEmrRdvDwSC2NkjNVWWLwnOAEyLiPSmlOyucryRJrbVaA8IgtjZIzVF116OoeH6SJLWO4WA8hgZpsTxHQZKkOTEgTM/QIM3fwoJCRARwNPDtlNJ1iyqHJEl1MhxUb8NDdjcsSHNQR1BIY463FXA2sLGmckiStBCGg/p16tjAINWnjh30b0XEd4DvAN8FvptSurrPeD9XPG+soQySJC2EIWG+bF2Q6lNHUHgQcCDwrM4bEbEeuLh4XAHcCjyjGHxNDWWQJGnuDAmLYeuCVI86gsJp5LDwUOB+xXs7AIcBh/YZ/yM1lEGSpLkxIDSDrQtSteoICn/QuY9CROxMDgydx4OBfYDtgeuAjwF/WEMZJEmaC0NCs9i6IFWn1pOIU0o3AucVD0mSWsWQ0FwGBml2my26AJIkLZsND9ndkLAk/J6k6RkUJEmagDuey8dgJ02nyqDw+8BZeLlTSVJLubO53Pz+pMlUdo5CSumvq5qXJElN4g5me3jugjQ+ux5JkjSEIaGd/F6l0QwKkiT1Yb/29vP7lYar9fKokiQto7bsQN54wFa1zXvny+6sbd7z5E3apMEMCpIkdVnWkFBnKJhkecsYIAwLUn8GBUmSCssSEuYdCiaxrAHCsCCVGRQkSaLZIaHJwWBc3Z+hqaHBKyJJmzIoSJJWvSaGhDaEg0GaHhpsXZAyr3okSVrVmhYSbjxgq1aHhF5N/bxNWy+kRbBFQZK0ajVpZ7CJO8vz1MRWBlsWtNrZoiBJWpWaEhKaekR9kZpUJ01ZT6RFMChIkladJuz8NWlnuKmaUkdNWF+kRbDrkSRpVVnkTl8TdnqXURO6JdkNSauRLQqSpFVjUSGhKUfG22CR9bjhIbvbuqBVxaAgSVoVFrGDZ0Cox6Lr1bCg1cKgIElqvXnv2C16R3a1MCxI9TIoSJJaaxFdRdoUENbvu/Fnj6YyLEj18WRmSVIrLaqrUVNUvXPfO781VzbnWOONB2zVmHsvSG1iUJAktc5qDQnzPPLftOCwqLDg1ZDUZgYFSVKrrLaQ0JRuQU0IDp3vYd6BwbCgtjIoSJJaY7WEhKaEg2H6lXFe4WERrQuGBbVRczoYSpK0ZOYZEpbhxOJR5ln+RQQ4T25W29iiIElqhXnupM1rJ3SZQ8Ew6/fdOJfWBU9ylmZji4Ikaem1LSQse8vBOOb1Ged9TwtbFdQmBgVJ0lJrY0hYTeYZGObFsKC2aF1QiIg9IuLMiLgmIu6IiHURcXpE7DzDPB8XEfdERIqIP62yvJKk6bUpJKyGVoRh5vH5DQvSZFoVFCJiX+Ai4IXA14DTgB8ArwK+GhG7TjHPNcD7gNsqLKokaYnMIyQoqzswGBak8bUqKADvAHYDXplSempK6bUppaPIgeEA4E1TzPOtwI7An1dXTEnSrOa1E1bnjuVqb0UYpi1hQVpmrQkKEbEPcAywDnh7z+A3ALcCJ0bEdhPM8ynk1olXAtdUU1JJ0rKoOyRouDqD1LzCgq0KWmatCQrAUcXz2SmlTbYqKaX1wPnAtsAh48wsInYD3g2clVL6QJUFlSTNZh47X4aE5rC+pMVo030UDiieLx8w/Apyi8P+wBfHmN+7yEHqpdMWKCIuGjDowGnnKUlaXovY4d1+75sqnd8tV+1Y6fzG1am7Ku+/MK/7LHjXZi2rNgWFzpZr0Bax8/5Oo2YUEScDTwF+LaX0PxWUTZJUkWVtTZhnSKg6HAyb97yDQ9U3a/OmbNJgbQoKo0TxnIaOFLEXcDrw0ZTSP82ywJTSwQOWcRHwiFnmLUmqxzKHhDoDwjjLnFdomNednatkq4KWUZuCQmdLNWgrtUPPeIOcCWwAXlZFoSRJ1am7NaHqkDCPgLCIcDDIPENDlWHBVgWpv+WK48NdVjzvP2D4fsXzoHMYOh5BvsTqdcUN1lJEJOC9xfDXF++dNVtxJUltVndI2H7vmxoVEnp1yldnGaus43lcBckrIGnZtKlF4Zzi+ZiI2Kz7ykfFTdMOJ7cUXDhiPu8nXx2p137A44BvkW/q9s2ZSyxJGtsytSbUGRKaHA4GWUT3pGnYsiBtqjVBIaV0ZUScTb6y0cuBt3UNPhXYDnhnSunWzpsRcWAx7aVd83llv/lHxEnkoPCZlNIfVv4BJEkL0/SQsIzhYJDt976p0rCwbOcreK6ClklrgkLhZcAFwBkRcTRwCfAY4Ehyl6PX94x/SfEcSJIaa1m6bFQdEtoUELo1OSzYqiCtWJ4IPoaU0pXAI4G15IDwGmBf4Azg0JTS9YsrnSSpiapqTagyJDT9/IMqVP0Zl+l8hWUJvlLbWhRIKf0IeOGY447dkpBSWksOIJKkOapzp2oeJ7BOqu0BoVeVrQvL1g1Jajp/TZIkzaiqo9mrLSR0NLFlwVYFyaAgSWqwZWhNMCRUo8quSMsSFqSmMyhIklYdQ8L4jtnzMo7Z87LRI1akaWGhTrYqqOlad46CJKkdmr4TtSwhoaqd/O75nP3DAyqZ5yBVnbdQxTkLXgVJq5ktCpKkVaWK1oRlCAl1tgTMo5WhSVd+qrMLUtMDsVY3WxQkSY1T185TU/qc17UDPM8uQr3Lq6uVYdbWBa+EJE3PoCBJ0gRmbU2oIyTMOyCMKkPVoaEJYaHOLkjerVlNZcSWJDVKk1sTmhYS5n2i8bjqKNesdVdFd7GmtEhJ82KLgiSp9doWEpoYDvqpupWhU4fTti40uWXBVgU1kS0KkiTVbDWGhF7LWm5pNTMoSJIao6lXgGnKNfmXfWe7CeVvyncpLQODgiRJWipNuWyq1HYGBUlSqy36BNSqdmqbcDS+Cm35HNJqYFCQJKnh3LlulrrCZ1O73mn1MihIkjTELH3aq2hNaGNIqOIzzVK3nqcgjcegIEmSJKnEoCBJUkO1sTWho82fTWoLg4IkSTWYtduRO9KjefUjqV4GBUlSI9RxIueir3ik5pr1PAXXLa0GBgVJkgbwpNd62WpS5pWP1CQGBUmSKma3o/mx+5FUH4OCJElaGEOR1FwGBUmSGsQd5/mxa5k0nEFBktRKs55sOu1OpF1hJjdrOFpUnXtCs9rOoCBJWjhP4Gy2E3e+YNFFkLQAWyy6AJIkKVtkt6NRYeDEnS/gH248rLblH7PnZZz9wwNqm/8y2fCQ3dnm4h8vuhiSQUGSpNVkltaBusPCIqzfdyNrrrSDhdSPQUGSpB6LOD+hrtaEqrsN1RkWZmlV2H7vm7jlqh0rLpG0uhmhJUlqqbrOLfCchRWe0Kw2MyhIklpn2Xbe6mhNqHtn3rAgtZ9BQZK0UG254lGTLos6r534OpYzS2ia9jvwfgpSfwYFSZK6uNM4GVsW6tGWAK3lZlCQJKlF2rDj7t2ppWYwKEiS1GXZL5U578uX1rE876eA91FQIyz31lCStPTq2CHa+bI7K5/nKLNcmrPqHeN/uPGw2gNDXcuYpS6m/Q6WPRxKdfE+CpIk9Vhz5WZzP1fh7B8eUHmXm94d+Sq6JdUZQJaxJWERoVSaF4OCJEkVueWqHRt19aNek+7kdweLulsoZg0JbWpNsNuRmqJ5vw5J0qrTxO5H0+5ANqkL0qw63YvaGhJmZWuC2s6gIElSgzQtLLRZE1sTpCbxFyJJaq1lbFWA1RUWlrU1oS52O1KTGBQkSY3Qth2ktu3A1mGRgWjW1gS7HWk1MChIkjTEorqnnP3DA1rdslDFZzOMSfUyKEiSWm2RR36r2JFtY1gwJPTXtlY1LT+DgiRJI8zSqlBVWGhLYGjC57DbkTQeg4IkqTHaekS1qqPfTdjJnkVV5Z+lPr3SkTQ+fy2SpNar4ghwU3Ywl7V1YRnL3E9drQltDclabs3Y6kmS1HJV96lflh3vqoONrQnS/PiLkSStCk1oVagjLDQ1MNRRtjaewCw1mUFBktQobe+CUcfObpPCQl3hZdZ6q6I1wW5HWm22WHQBJElaJmuu3Iz1+25cdDFKOjvnx+x52UKXXwdbEqTFsEVBkrRqNOWylnXu+M67daHu7k9V1JXnJkjT8ZcjSWqcpnfFqGLHc9nDwjzOj2hSS4LdjrQa2fVIkqQFueWqHdl+75tqmXf3TnyV3ZGadD7EOGxNkKZnUJAkrSo7X3YnNx6w1czzqepchTrDQse0O/fdAWPeAWE1tCZITWdQkCRJfS2q9aCqkND01gS7Hanpmv0LkiStWnXuRFV1hLiqHdEmHT1ftNUSEqRlYIuCJEkzqLILElB7N6SmqjIsVRkSPIlZq5lxW5LUWHW3KjStZQHyDvNqamGo+vNW9V1UuX70MiRoWRgUJEmNVvdOVZVhwcAwvjoCQpUhoS6GBC0Tg4IkqfGWJSxA9X3jOzvUbQkOdXyOZehqBIYELR/PUZAkLYVtLv4xGx6ye23zr+qyqbCy41rFuQu9eneym3xOQ93BpupQZkiQNmVQkCQtjWUKC1BvYOjo3hlfdGiYZ4uHIUGqn0FBkrRUli0swHwCA8y3tWFR3aDquOypIUHqz6AgSVo68wgLwNIGho42nNPQUdd9EbzrsjSYQUGStJTqDgtQT+sCzD8wLKs6b5o2j4Bga4KWnVc9kiQtrXnsiNV5Pf2qL6naFnXXiyFBGo8tCpKkpTaPlgWorzsS2MIA9bYedMyrm5EhQW3hYQxJ0tKb545ZnTubq7GFYR6fuc5WoV6GBLWJLQqSJE2oztYFKB9db1tLwzzD0DxPVjYkqG0MCpKkVphXF6RudZ3s3KvfjvWyhYd5t5TM+2pGhgS1kUFBktQanZ21eQaGulsXBmlqq8Oiu04t4nKnhgS1lUFBktQ6i2pdgPkHho55tjosOgz0s6j7IRgS1GYGBUlSKy2idQHm1x1pHE3coa/aIm+YZkhQ2xkUJEmttsjWBVhcC0PbGRCk+hkUJEmtt6jWBTA0VG2RAQEMCVpdDAqSpFVjEa0L3QwNk1t0MOgwIGg1MihIklaVRbYudOvdATY4rGhKOOgwJGi1al1QiIg9gDcCxwK7Aj8BzgJOTSndOMb02wFPBX4FeARwf2AjcBnwIeBtKaVmbcEkSRNbdOtCr9Xc2tC0YNBhQNBq16qgEBH7AhcAuwGfAi4FHg28Cjg2Ig5PKV0/YjaPBT4A3ACcQw4ZuwBPAv4GOD4ijk4p3V7Pp5AkzUtTWhd6tb21oanBoJshQWpZUADeQQ4Jr0wpva3zZkS8Bfht4E3AS0fM47+B5wEf7W45iIg1wLnAYcDLgTdXWnJJ0sI0rXWhV78d62UKD8sQDDoMCNKKSCktugyViIh9gCuBdcC+KaWNXcPWkLsgBbBbSunWKZfxHOAfgU+nlJ40Q1kv2mHL+zzisN1OmHYWkqSaNDkwjLLI8LBMYaAfA4La6oJrP8zNd133Hymlgyedtk0tCkcVz2d3hwSAlNL6iDgfOAY4BPjilMu4q3i+e8rpJUkN19TuSOMYZ2d91jCx7IGglwFBGqxNQeGA4vnyAcOvIAeF/Zk+KJxcPH9+nJEj4qIBgw6ccvmSpDlZ5sAwTNt29GdhSJCGa1NQ2LF4vmnA8M77O00z84j4LfKVlL4FnDnNPCRJy6etgWE1MyBI42lTUBgliueJT8qIiOOB08knOj89pXTXiEnyggb0BStaGh4xaTkkSYvTvXNpaFg+hgNpcm0KCp0Wgx0HDN+hZ7yxRMRTgQ8D1wJHppR+MF3xJEltYSvD8jAgSNNrU1C4rHjef8Dw/YrnQecwlETEM4EPklsSjkopXTF98SRJbWNgaCbDgVSNNgWFc4rnYyJisz6XRz0c2ABcOM7Mikuhvh/4MbYkSJKGMDA0gwFBqtZmiy5AVVJKVwJnA3uRb4jW7VRgO+D93fdQiIgDI6J0BaKIeAHwD8APgccZEiRJ49jm4h+7szpnnTq33qXqtalFAeBlwAXAGRFxNHAJ8BjgSHKXo9f3jH9J8dw50ZmIOJJ8VaPNyK0UL4yInsn4aUrp9MpLL0lqhd6dVlsaqmUokOajVUEhpXRlRDwSeCP5UqZPJN+R+Qzg1JTSDWPM5gGstLScPGCcq8lXQZIkaaR+O7aGh/EZDKTFaFVQAEgp/Qh44ZjjlpoKUkprgbXVlkqSpE3Z6tCfoUBqjtYFBUmSltFqDA6GAqnZDAqSJDXQODvRTQ8TBgFpuRkUJElaUnXtiA8KIO74S6uLQUGSJG3CQCAJWnQfBUmSJEnVMShIkiRJKjEoSJIkSSoxKEiSJEkqMShIkiRJKjEoSJIkSSoxKEiSJEkqMShIkiRJKjEoSJIkSSoxKEiSJEkqMShIkiRJKjEoSJIkSSoxKEiSJEkqMShIkiRJKjEoSJIkSSoxKEiSJEkqMShIkiRJKjEoSJIkSSoxKEiSJEkqMShIkiRJKjEoSJIkSSoxKEiSJEkqMShIkiRJKjEoSJIkSSoxKEiSJEkqMShIkiRJKjEoSJIkSSoxKEiSJEkqMShIkiRJKjEoSJIkSSoxKEiSJEkqMShIkiRJKjEoSJIkSSoxKEiSJEkqMShIkiRJKjEoSJIkSSoxKEiSJEkqMShIkiRJKjEoSJIkSSoxKEiSJEkqMShIkiRJKjEoSJIkSSoxKEiSJEkqMShIkiRJKjEoSJIkSSoxKEiSJEkqMShIkiRJKjEoSJIkSSoxKEiSJEkqMShIkiRJKjEoSJIkSSoxKEiSJEkqMShIkiRJKjEoSJIkSSoxKEiSJEkqMShIkiRJKjEoSJIkSSoxKEiSJEkqMShIkiRJKjEoSJIkSSoxKEiSJEkqMShIkiRJKjEoSJIkSSoxKEiSJEkqMShIkiRJKjEoSJIkSSoxKEiSJEkqMShIkiRJKjEoSJIkSSoxKEiSJEkqMShIkiRJKjEoSJIkSSoxKEiSJEkqiZTSosuw6kTE9ZuxxS7bb7nzoosiSZKkFrvlrhvZyN03pJR2nXRag8ICRMQdwObAtxddlpY5sHi+dKGlaB/rtR7Waz2s13pYr9WzTuthvZbtBdycUtp70gm3qL4sGsPFACmlgxddkDaJiIvAeq2a9VoP67Ue1ms9rNfqWaf1sF6r5TkKkiRJkkoMCpIkSZJKDAqSJEmSSgwKkiRJkkoMCpIkSZJKvDyqJEmSpBJbFCRJkiSVGBQkSZIklRgUJEmSJJUYFCRJkiSVGBQkSZIklRgUJEmSJJUYFCRJkiSVGBQqEBFbRsSrIuK9EfGtiLgzIlJEvGiGeR4WEZ+NiBsi4raI+E5EvDoiNh8yza9GxLkRcVNE3BIR/x4RL5i2DE0xTV30mccpxXcy7HFlzzRHjBj/L6r/tPNRRZ0W8xlWPxcOmc51dfA8do+IV0TE5yJiXUTcERHXR8S/RsTxA6ZZ+nU1IvaIiDMj4priM6+LiNMjYucJ57NLMV2n7q4p5rtH3ctuolk/W0RsFxHPjYgPRsSlEXFrRKyPiG9ExGsiYqsB0021bVgGVawvxfZvWB1tPWC6X4iIf4qIayPi9oi4LCJOjYhtqvuEi1HBujpqO9h53L9nutauq1XYYtEFaIntgNOL1/8D/Ddw/8GjDxcRTwE+DtwOfAS4AXgScBpwOPDMPtP8FvA24HrgA8CdwDOAtRFxUErpd6ctzyJNUxcDnDtk2JOARwCfGzD8vAHTf2XMZTdKhXXacTWwts/7/zVg+a6rw70C+H3gKuAc8vbkAcDxwC9FxGkppd8ZMO1SrqsRsS9wAbAb8CngUuDRwKuAYyPi8JTS9WPMZ9diPvsDXwI+DBwIvBD4lYg4NKX0gzqW3UQVfbbHkn+nN5DXx7OAXcjr9t8Ax0fE0Sml2/tMO9G2YRnUsL6cOuD9u/ss+zHk9XpL4GPAj4CjgD8Gji6+hzsmWHZjVFSv6xhcnweRt6HfSyn9qM/w1q2rlUkp+ZjxAWwFHAf8fPH3KUACXjTFvHYArgXuAB7Z9f7W5B9RAk7omWYv8s7J9cBeXe/vDHy/mObQRdfTPOpiimVsTt7YJuChPcOOKN4/ZdF10dQ6LcY/d4LxXVdHz+t44PF93n8QcFMxr4N7hi31ugr8S1H+V/S8/5bi/b8bcz7vLMZ/S8/7ryze/3xdy27io4rPBjwceC6wVc/7a4CLivm8ps90E20bluVR4bp6LpAmWO7mwH8Wy3hy1/ubkUNDAl676PpZdL0Omf+Hivm8ss+wVq6rlX03iy5AGx/MFhROLqZ9X59hRxXDzut5/43F+6dOMr+mP6apiymW8aRiPl/tM2ypd77mUaeTbmBdV2dezrv67Zgt87oK7FOU/Spgs55ha4BbgFuB7UbMZzvgtmL8NT3DNivmn4B9ql52Ex/z+GzAc4pl/HOfYa3b+aqyTpk8KAzcjnSVax0Qi66nRdbrgPnvSj5AdRuwc5/hrVtXq3x4jkLzHFU8f77PsH8jr+iHRcS9xpzmcz3jLJNp6mJSv1E8v2vIOA+MiN+KiNdFxMkRsd8My1u0Oup0p6JeXhcRL4+IQ6ZcvuvqaHcVz6VuCYVlXFc7dXd2Smlj94CU0nrgfGBbYNh6BXAosA1wfjFd93w2AmcXfx5Zw7KbaB6fbdT6OMm2YRlUXqcR8WsR8dqI+J2IOG7INmLgNibl7nSXk7so7jPushuk7nX1JOBewEdTSjcOGKdt62plPEeheQ4oni/vHZBSujsirgIeTN4YXDLGND+JiFuBPSJi25TSbTWUuS7T1MXYImJ3cpexm8h9ygd5bvHonvbjwIuHbHSaqo46fRjw991vRMS3gRNTSt+dYPmuq0NExA7A08lHv84eMNoyrqsD665wBXAM+byDL844H4r5VL3sJprHZzu5eO4XkGGybcMyqKNOP9zz97UR8fKU0semWPb+xePKAeM0Vd3raufCMu8cMk7b1tXK2KLQPDsWzzcNGN55f6cpptlxwPCmmqYuJvEicr/PDwzYKb0OeC35JKg1wH3IweKb5B22f46IZfsNVV2nbyGfqHsfch09itxf9mHAl4owNs3yXVe7REQA7wHuC/xtSqk3bCzzulpV3dW57Zx2G7NIda+TvwUcC3wLOLPPKJNuG5ZBlXX6KXLX1z3ILWEHAn9eTPuRiDiuxmU3TW2fLSIeT67b76WULhgwWhvX1co09R/H3BWX4RrnslqdxwcWVdTiOdU8TSVqrtepP1ex09Q5Gta321FK6Xsppb9MKV2cUrolpfS/KaXPk/uDX0XesDxp0mXPqkl1mlJ6TUrpgqJubkkpfSOl9Ezy1X/uDUx6BSPX1f7eTL5q0peB0hWPmrquVqSqdWKptp1zMMv283jylf7+G3h6Sumu3nFq2DYsg7HrNKV0Wkrp0ymlH6eUbk8pXZZSeh3wGvK+2Z/VtewlNMtn63QvHtiasErX1bHZ9WjFleSTXcZ1TU3lGHVEdYee8Tqv711M0+/yYZ1pbp65dJObpV6nqYtxHQfsCVyYUvrOJBOmlG6OiA8CrwceRz4yNE9NrdNuf0c+kv24nvddVycUEX8N/Db5XIdfSRNc/rAB6+o4qqq7abedVSy7iWr5bBHxVHJ3mWuBI1PP5WbHMGjbsAzmsb68h3xJ5YdHxJqu821cVydfV3chr2sbgH+YolzLvK5WxqBQSCkdvegyFC4DHknui3dR94CI2ALYm3zi2A96prl3Mc1Xe6b5efLVQP5rEX2+Z6zXaepiXCOPMoxwXfG83ZTTT63BddptUP24rk4gIk4DXk2+fv2vTlkvC1tXx3RZ8bz/gOGdE7IH9V+eZT5VLbuJKv9sEfFM4IPkloSjUkpXjJikn6avj8PUvr6klG6PiPXkS0ZvB3SCguvq5J/tBeSTmN+XUvrpFOVa5nW1MnY9ap4vFc/H9hn2OPKZ/xf0HFUcNs1xPeMsk2nqYqSIuB/wkJmLlwAABnRJREFUK+SjE/80Zdk6V0SYdYd63mqp0z4G1Y/r6hgiezs5JPwruSVh2vDU9HX1nOL5mN7zKCJiDbnb1AZg1B1SLyzGO7yYrns+m5FPhuxeXpXLbqJKP1tEPId8LfpryPf5mCYkQPPXx2FqX18i4gBySFgP/G/XoIHbmIjYh7yTfTXWa7cXF8/Drmo4zDKvq9VZ9PVZ2/hgjPsokJvYDqS4SVvX+zuQU+wkN1zbm/bexGrSuti2qNc9h8z3j4pp3zZi+YfTc03n4v3nARuLcu017udpwqPKOiXfzbp0XWvgoeR/cAl4juvqxPUawLuLaT4LbD3G8pd6XWXCmy0V9XZgn/l0brj25p73veHabPX6AuAe8g7TA8ZY7sTbhmV5VFGn5Kuf7d5n3vfu2l68q2fYsBuufZRVdsO1Qetq1/DHFtN9d8RyW7uuVvWIokI0o4h4LXnFhXwny4eRf/Cdoy5fSSm9p2v8k4D3kpvETuqZ11PJZ9zfTu4LegPwZPIlxD4GPCv1fHER8QrgDPIO2EeAO4FnkK+o8OaU0lKejDNpXUTEEeSjE+ellI7oM7/NKP7Zke/EPPCyZxGxjrwRvoB8G/etyVdDeDS5G8mLU0prZ/yIc1dVnUbEWvJdhL9Evrv1HeTfwLHkf2rvBl7iujpxvb6BfLBhA/mE0Tv7LO5bKaWzuqZZxxKvqxGxL7nsu5HPo7gEeAz5ngeXA4ellK7vGj/fJSml6JnPrsV89ievl18j39H6KeQ+9YellK7smWaiZS+TKuo1Io4EvkBev84k/9Z7/TSldHrXNGuZYtuwDCqq05PI5yKcRz436gbyOXNPJB9E/Abwy6mnu0xEPIZcp1uStyk/BI4md3s8Hzg6zd4avBBVbQO6hv8D+UDJK1NKbxuy3LW0dF2tzKKTSlseFHdZHPJY2zP+Sf3e7xp+OPlo4o3kHYbvkk9o3HxIGZ5E3vCsJ9/F8OvACxZdNxXU7dh1wcodas8dMK/jGHAn5j7j/j6528ePiuXeTt6ovxd42KLrZdF1CjwV+AS5JeBm8g7tT4B/puuIl+vqxPW6dsS2pN/2ZOnXVeD+RXl/UqxLVwNvBXbpM25iwF1tgV2K6a7uWifPBPaoYtnL9pi1Xln5XzXssa5nmqm3DcvwqKBODyp+598lHzC5ixwWvgy8AthqyLJ/gdyC8L/kndrLgVOBbRZdL4uu165hOxfbwduAnUYss9XrahUPWxQkSZIklXgysyRJkqQSg4IkSZKkEoOCJEmSpBKDgiRJkv7/9u5YRaorjuP4768BsVAwbUi1aCGiJiCBVOnzANpIEPIipgv4BLJNAkmeIeQVEkwwgultBAvRIgrqSTH3kiFnuZFlZ+ew8/k0d+7MKf7tl3PuXOgIBQAAoCMUAACAjlAAAAA6QgEAAOgIBQAAoCMUAACAjlAAAAA6QgEAAOgIBQAAoCMUAACAjlAAYEhVtV9Vrap+WVhzrqreTeu+PM75AE46oQDAqK5O198X1lxPUtPnPzY7DsBuEQoADKeqTiW5Mt0+WFh6fbo+a6092exUALtFKAAwoktJzk6f/29HIbGbAHDkhAIAI5qPHf2d5K+FdZ9M16WYAOAQhAIAI7o2XR+21t4etKCqPkhyeboVCgBHTCgAMKI5FJYC4HKSM++xDoBDEAoAjGg+erT0IPN87OhVksebHQdg9wgFAIZSVReSfDzdLu0UzKHwqLX2ZrNTAeweoQDAaK6tfX64sO6L6erYEcAGCAUARjMfO3qb1bGiTlXdyPs9xwDAIQkFAEYzB8DpJJ/998eqOp/k/tpX3qEAsAFCAYDRrB89+qGqblbVXlVdrKo7SX7Lvy9aS5K9qvr8eEcEOPmqtbbtGQAgSVJVp5O8zOqtzPtJvj5g2Ysk95J8s/bd9621rzY/IcDusKMAwEguZRUJSXI3ye2sHmh+neRpkh+TfJrk2yQ/Z/Xm5l+TfHfskwKccHYUABhGVd1K8lOS5621C9ueB2CX2VEAYCTzPx492uoUAAgFAIYyP8j851anAEAoADCUORTsKABsmVAAYAhV9WGSj6ZbOwoAWyYUABjF+vsThALAlvnXIwAAoGNHAQAA6AgFAACgIxQAAICOUAAAADpCAQAA6AgFAACgIxQAAICOUAAAADpCAQAA6AgFAACgIxQAAICOUAAAADpCAQAA6AgFAACgIxQAAIDOP8FbbGg7kd0wAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 287,
       "width": 389
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "Xfact = np.mgrid[-1:1:0.05, 0.01:1:0.05]\n",
    "Zfact = np.apply_along_axis(q, 0, Xfact, D=D, E_mu=E_mu, E_tau=E_tau)\n",
    "\n",
    "title = r\"Factor posterior distribution $q^*_\\mu(\\mu)q^*_\\tau(\\tau)$\"\n",
    "plt.title(title, fontsize=15)\n",
    "plt.contourf(*Xfact, Zfact)\n",
    "plt.ylabel(r\"$\\tau$\", fontsize=13);\n",
    "plt.xlabel(r\"$\\mu$\", fontsize=13);"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
